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' We consider one dimensional interacting bose-fermi mixture with equal masses of bosons and 

fermions, and with equal and repulsive interactions between bose-fermi and bose-bose particles. 
Such a system can be realized in current experiments with ultracold bose-fermi mixtures. We 
OA ' apply the Bethe-ansatz technique to find the exact ground state energy at zero temperature for any 

value of interaction strength and density ratio between bosons and fermions. We use it to prove 
q . the absence of the demixing, contrary to prediction of a mean field approximation. Combining 

exact solution with local density approximation (LDA) in a harmonic trap, we calculate the density 
profiles and frequencies of collective modes in various limits. In the strongly interacting regime, we 
00 , predict the appearance of low-lying collective oscillations which correspond to the counterflow of 

the two species. In the strongly interacting regime we use exact wavefunction to calculate the single 
particle correlation functions for bosons and fermions at low temperatures under periodic boundary 
conditions. Fourier transform of the correlation function is a momentum distribution, which can be 
■ measured in time-of-flight experiments or using Bragg scattering. We derive an analytical formula, 

which allows to calculate correlation functions at all distances numerically for a polynomial time in 
the system size. We investigate numerically two strong singularities of the momentum distribution 
for fermions at kf and fe/ + 2fei,. We show, that in strongly interacting regime correlation functions 
change dramatically as temperature changes from to a small temperature ~ Ef/f <C Ef, where 
Ef — (irhn) 2 1 (2m) , n is the total density and 7 = mg/(h 2 n) ^> 1 is the Lieb-Liniger parameter. 
A strong change of the momentum distribution in a small range of temperatures can be used to 
perform a thermometry at very small temperatures. 
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q ■ I. INTRODUCTION 

^v^j I Recent developments in cooling and trapping of cold atoms open exciting opportunities for experimental studies 
J> . of interacting systems under well controlled conditions. Current experiments^ can deal not only with single com- 
t-H ' ponent gases, but with various atomic mixtures. Using Feshbach2ii resonances and and/or optical lattices^ one can 
tune different parameters, and drive the systems towards strongly correlated regime. The effect of correlations is 
most prominent for low dimensional systems, and recent experimental realization 7 ^ of a strongly interacting Tonks- 
Girardeau (TG) gas of bosons opens new perspectives in experimental studies of strongly interacting systems in IDS. 
In this article we investigate bose-fermi mixtures in ID, using exact techniques of the Bethe ansatz. Some of the 
I results presented here have been reported earlier—. 

Most of the theoretical research on bose-fermi mixturesii so far has been concentrated on higher dimensional 
systems, and only recently ID systems started attracting attention. Several properties of such systems have been 
investigated so far, including phase separationi2iiii4, fermion pairing^, possibility of charge density wave (CDW) 
1 1 formation 16 and long distance behavior of correlation functions^ 7 -. 

A ID interacting bose-fermi mixture is described by the Hamiltonian 

« , 

O ■ r L fi 2 

O ■ H = 



/ dx{— d x Vld x V b + dMd x ^ } )+ dx(-g bb ^l^ b ^ b + g bf ^M^ f ^ b ). (1) 

J Zm b Zrrif 1 J I 1 

Here, ty b ,tyf are boson and fermion operators, m b ,rrif are the masses, and g bb ,g b f are bose-bose and bose-fermi 
interaction strengths. The model (JJJ is exactly solvable, wheniS 

m/ = m b = to, g bb = g bf = g > 0. (2) 

It corresponds to the situation when masses are the same, and bose-bose and bose-fermi interaction strengths are the 
same and positive. Although conditions (J5J are somewhat restrictive, the exactly solvable case is relevant to current 
experiments (the experimental situation will be analyzed in detail in section lVlII ) and can be used to check the validity 
of different approximate approaches. Model under conditions (0) has been considered in the literature beforei^, 
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but its properties have not been investigated in detail. After the appearance of our initial report, two additional 
articlesi2*i*i used Bethe ansatz to investigate the same model. We use the exact solution to calculate the ground state 
energy and investigate phase separation and collective modes at zero temperature. For strongly interacting regime, we 
calculate single particle correlation functions, and consider the effects of small temperature on correlation functions 
and density profiles. 

The article is organized as follows. In section^ we review the Bethe ansatz solution for bose-fermi mixture and 
compare it to the solution for fermi mixture. In section ITTT1 we obtain the energy numerically in the thermodynamic 
limit. We use it to prove the absence of the demixing under conditions @, contrary to prediction of a mean fielcU^ 
approximation. In section llVl we combine exact solution with local density approximation (LDA) in a harmonic trap, 
and calculate the density profiles and frequencies of collective modes in various limits. In the strongly interacting 
regime, we predict the appearance of low-lying collective oscillations which correspond to the counterflow of the two 
species. In section we use exact wavefunction in the strongly interacting regime to calculate the single particle 
correlation functions for bosons and fcrmions at zero temperature under periodic boundary conditions. We derive an 
analytical formula, which allows to calculate correlation functions at all distances numerically for a polynomial time 
in system size. In section IvTl we extend the results of sectionEJfor low temperatures. We also calculate the evolution 
of the zero temperature density profile at small nonzero temperatures. We show, that in strongly interacting regime 
correlation functions change dramatically as temperature is raised from to a small value. Finally in section IVIII we 
analyze the experimental situation and make concluding remarks. 



II. BETHE ANSATZ SOLUTION 



In this section we will briefly review the solution 18 of the model Q under periodic boundary conditions and compare 
it to the solution of Yang of the spin-i interacting fermions^SI, for the sake of completeness. More details on Yang's 
solution can be found in 22 i 23 i 24 i 25 i 26 . 

In first quantization, hamiltonian can be written as 

N d 2 

H = -22 7rs + 2cJ2S(x i -x j ),c>0. (3) 

Here we have assumed m = 1/2 and h — I, to keep contact with the literature on the subject. Later in the discussion 
of the collective modes we will introduce the mass of atoms, but it should be clear from the context whether we have 
assumed m = 1/2 or not. c in is connected to parameters of via 

Wave function is supposed to be symmetric with respect to indices i = {1, M} (bosons) and antisymmetric with 
respect to i — {M + I, A} (fermions). On the first stage, Yang's solution doesn't impose any symmetry constraint 
on the wavefunction. On the second stage, periodic boundary conditions are resolved with the help of extra Bethe- 
ansatz. This idea has been generalized by Sutherland 2 ^ for the case of N-fermion species. The results presented here 
can be simply derived from Sutherland's work. 

In Yang's solution, one assumes the generalized coordinate Bethe wavefunction of the following form: for < xq 1 < 
xq 2 < ... < x Qn < L 

* = Y}Q,P¥^ XQ ^E = Y J kl (5) 

P i 

where k\, ...,/cjv is a set of unequal numbers, P is an arbitrary permutation from Sjv and [Q,-P] is A! x AM matrix. 
Let's denote the columns of this matrix as A! dimensional vector £p. Delta function potential in J2J is equivalent to 
the following boundary condition for the derivatives of the wavefunction: 

d d d d 

("^ — )^x j= x k +0 - "o — )^x j= x k -0 = 2c$> X]=Xkl (6) 

oxj OXk oxj oxk 

and the continuity condition reads 



*x 3 =x fc +0 = ,. ::• (7) 
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Suppose Q and Q 1 are two permutations, such that Qk = Q' k , for k ^ {i,i + 1}, and Qi = Q'i + i,Q'i = Qi+i- 
Similarly, P and P' are two permutations, such that P k — Pf., for k ^ {i,i + 1}, and Pj = P- +1 ,P' i = Pi+i- To 
satisfy © and J2J for xq ( = xq i+1 independently of other x, one has to impose two conditions for four coefficients 
[Q,P],[Q',P} 7 [Q,P'},[Q', Z P'}. Using these two conditions, we can express [Q , P'] , [Q' , P'} via [Q, P],[Q' , P]. These 
requirements can be simply written as a condition between £p and £p< : 
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M+l 



(8) 



y operators are defined as 



Y 



1 + A, 



1 + A, 



-P 



(9) 



where 



hi 



and P/ m is an operator acting on a vector £p which interchanges the elements with indices Qi and Q m . Using Y 
operators one can express any £p via £o? where is a column for P — identity. However, arbitrary permutation P 
can be represented as a combination of neighboring transpositions by different means. Independence of the final result 
on a particular choice of neighboring transpositions can be checked based on the following Yang-Baxter Relations: 



j,k i,k i,j 



■*jrct,b-*rct,b i 

yb.cya,b-yb,c 
i,j i,k j,k ' 



(10) 

(11) 



Operators Yp l ~^> exchange the momentum labels Pj and Pf+i, while P^. 
and Qi+\. It is convenient to define combined operator, which exchanges both labels: 



interchange relative position labels Qi 



P-Y 13 



i _ \ p 

1 + A;,- 



(12) 



Using this definition, periodic boundary conditions can be written as N matrix eigenvalue equations: 



Xj+l,jXj+2,j- 



..V v A , V, i 



-ikj L 



Co 



(13) 



The procedure outlined above reduces equations for AM x AM coefficients to Af eigenvalue equations for AM dimensional 
vector. Imposing some symmetry on £ simplifies the system further. If £ is antisymmetric with respect to particle 
permutations (fermions) , then P^ = — 1 and e lkjL — 1. The system of equations is the same as for noninteracting 
fermions, as expected. If £o is symmetric (bosons), Py = 1 and the system is equivalent to periodic boundary 
conditions of Lieb-Liniger modetf£. 

If one needs to consider two-species system, £o has the symmetry of the corresponding permutation group represen- 
tation (Young tableau). Instead of solving eq. I|13l) . it is convenient to consider the similar problem in the conjugate 
representation. If £o is antisymmetric with respect to both permutations of the first M indices and the rest A^ — M 
(two-species fermions), eigenstate in conjugate representation tp is symmetric with respect to first M indices and is 
also symmetric with respect to permutations of the rest N — M indices. Similarly, in conjugate representation for 
bose-fermi mixture with M bosons and N — M fermions tp should be chosen to be antisymmetric for permutations 
of M boson indices and symmetric with respect to permutations of A^ — M fermion indices. The periodic boundary 
conditions are (note the change of the sign in the definition of X^ compared to Xij): 



y' y' 



Xn,-)X[ 



.X'tp = V 



Ay Pij 



1 + Ai 



(14) 
(15) 



Since AM-dimensional vector tp has symmetry constraints, it has Cjjf inequivalent components, characterized by the 
positions yi of M spin-down fermions (or M bosons respectively). One can think of the components of the vector tp 
as of the values of the spin wavefunction, defined on an auxiliary one-dimensional lattice of size N. C^Jf independent 
values of tp correspond to CfJ values of the wavefunction of M "particles" with coordinates yi, living on this auxiliary 
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lattice (since <p is symmetric for N — M fermion indices, these are considered to be vacancies). Wavefunction should 
be symmetric with respect to exchange of two "particles" for two-species fermions, and antisymmetric for the case 
of bose-fermi mixture. To preserve the terminology of the two-species fermion solution for the case of bose-fermi 
mixture, later in the text we will always refer to the wavefunction on an auxiliary lattice as to "spin" wavefunction, 
although it has a direct meaning only for two-species fermion case. 

First, one can solve the problem for M — l 30 . In this case there is no difference between two-species fermions 
or bose-fermi mixture. It can be shown (detailed derivations are available in the appendix oS&), that in this case 
wavefunction in conjugate representation is 

where new spectral parameter A satisfies the following equation: 

N 



h - A + ic/2 

«=i 

Periodic boundary conditions simplify to 

al _ kj-A + ic/2 



n**" ^ 1 =1. (17) 

life -A- ic/2 1 ' 



e 



kj - A - ic/2 



j ={!,... ,N} (18) 



In an auxiliary lattice the wavefunction of one spin-deviate (or boson) F(A, y) plays the role similar to one-particle 
basis function e tkx of the original coordinate Bethe ansatz, spectral parameter A being the analog of the momentum 
k. 

In the case when M > 1, Yang suggested that the solution of eqs. I|14|l - H15|) again has the form of Bethe ansatz in 
the "spin" subspace: for 1 < y\ < yi < ... < yu < N 

M 

^ = Y / MR)l[F(A Rt ,y t ), (19) 

R z=l 

where Ai, Am is a set of unequal numbers, R is an arbitrary permutation from Sm- It can be shownSSaSi, that this 
ansatz solves I|14f) - (I15II for two-species fermion system, if 

MBf) = A Rt+1 - A Ri - ic 

A{R) A Ri+1 -A Ri +ic' 1 ) 

similarly to bosonic relations of Lieb-Liniger models. Here R and R' are two permutations from Sm such that 
Rk = R' k , for k 7^ {i, i + 1}, and Ri = R' i+ i, R[ — Ri+i- The set of A, k has to satisfy the following set of equations: 

n Y k l -A a +ic/2 j£ A - A a + ic 

For the bose-fermi mixture, if has to be antisymmetric for permutations of yi variables. This problem has actually 
been solved by Sutherland 27 , although he was interested not in bose-fermi mixture, but fermion model with several 
species. He has shown, that if one doesn't specify the symmetry of ip for y t variables and applies the generalized 
ansatz 

M 

<p = Yfe,R]l[F(A at ,yG i ) (23) 

R i=l 

for 1 < ya t < yc 2 < ■•• < Vg m — then columns of Ml x Ml dimensional matrix [G, R] are related similarly to JHJ: 

& = Y '£tlM- ( 24 ) 
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Y' operators are defined as 

yll,m K ij Plm ^ (25) 

For two-species fermions in conjugate representation P\ m — 1, and it is equivalent to (|2U[I . while for bose-fermi mixture 
in conjugate representation Pi m — — 1, and the answer is much more simple: 

Yff 1 = -1. (26) 

Therefore, "spin" part of wavefunction is constructed by total antisymmetrization of single "spin" wavefunctions, 
similar to Slater determinant for fermionic particles: 

tp = det(F(A i ,y j )). (27) 

Periodic boundary conditions for bose-fermi mixture are: 

n =1 -" ={1 -- M} - ^ 

i— 1 ' 

t tt fc-i ~~ Afl + ic/2 - , 

^ L = U k _ A -^ = < 1 - jV >- ( 2Q ) 

/3=1 3 p ' 

One can prove that all solutions of (|28|l - 1)29(1 are always real, which is a major simplification for the analysis of both 
ground and excited states (see Appendix . 
If one introduces function 

Q(k) = -2tan~ 1 (fc/c), (30) 

the system (|28|) - (|29f) can De rewritten as 

M 

fcj-L = 2ttI 3 + ^ 6»(2fc 3 - 2A (3 ), (31) 

AT 

27r/ Q = ^ 6>(2A Q - 2%). (32) 

3=1 

7j and I a are integer or half integer quantum numbers (depending on the parity of M and N), which characterize the 
state. The ground state corresponds to 

I a = (M - l)/2, —{M - 3)/2, (M - l)/2}, (33) 
Ij = {-(N - l)/2, —(N - 3)/2, (N - l)/2}. (34) 

In the thermodynamic limit, one has to send M,N,L to infinity proportionally. If one introduces density of k roots 
p(k) and density of A roots cr(A), <|28|) - (I29[1 simplifies to two coupled integral equations 

2np(k) = l+ [ / Ca ^ d \ 2 , (35) 



B c 2 + 4(A - fc) 2 

Acp(uj)duj 
QC 2+4(A-c)2 



MA)=/ J°± )(LJ ., 2 - 06) 



Normalization conditions and energy are given by 



N/L= / p(k)dk, (37) 
J-Q 

M/L= [ a(A)dA, (38) 

J-B 

rQ 

E/L= I k 2 p{k)dk. (39) 
J-Q 

These equations can be solved numerically and the results will be presented in the next section. Numerical solution 
of these equations allows to investigate the possibility of phase separation, predicted ini^. Combined with local density 
approximation, it can be used to investigate density profiles and collective oscillation modes in the external fields. 
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III. NUMERICAL SOLUTION AND ANALYSIS OF INSTABILITIES 



In this section we will solve the system of equations <|35f) - (|39|) numerically, and obtain the ground state energy 
as a function of interaction strength and densities. This solution will be used to analyze the instability towards 
demixingi2iiii4. 

Substituting l|36|) into (|35(l . and performing analytically integration over A, one obtains an integral equation for 
function p{k). Similarly to2£, it is convenient to redefine the variables before solving this equation numerically. Let's 
introduce the following variables A, x, y, 6 and a function 17(0;) according to 

c = \Q,uj = xQ,k = yQ, B = bQ, p{Qx) = g{x). (40) 

In new variables, integral equation depends on two parameters b and A : 

„ f 1 2Xg(x)dx . ,2(6-1/) ,2(6-x) 

tan - w+»> + tm - w+j) + * , oe * ± f - % 1 1 f +»>;,. (41) 

A A 2(x-y) A 2 + 4(6 - y) 2 A 2 + 4(6 + x) 2 K 



In new variables, (|37|I - H39|) become 



cL A , , . 

7 = ^= f i : w ■ ( 42 ) 



M _ /^(tan- 1 + tan- 1 ^1 )g( x )dx 

N it J g(x)dx 

N 3 N 3 J\x 2 g(x)dx 



(43) 



g = _ e( A,6) = ^ ;-; (44) 

J_ lff (x)da;J 



Integral equation (|41ll can be solved numerically as a function of two parameters 6 and A, applying Simpson rule for an 
integral approximation on a grid Xi = — 1 + (i— l)/n, i = {1, 2n + 1}. This gives a system of 2n + 1 linear equations 
for discrete values g{xi), which can be solved by standard methods. Using ^42 (1 -1)44 (1 . one can obtain parametrically 
three functions 7(A, 6), M/N = a(A, 6), e(A, 6). After that one can numerically inverse two of them A(7, a) and 6(7, a), 
and obtain function e(7, a). Resulting function is shown in fig. ^ When a = 0, system is purely fermionic, and 
noninteracting. When a = 1, the system is purely bosonic, and numerically obtained energy coincides with the result 
o^£. If 7 = 0, bosons and fermions don't interact, and e(-f, a) — (7r 2 /3)(l — a) 3 . 

An interesting case, where one can analytically find the dependence of energies on relative densities is Tonks- 
Girardeau (TG) regime of strong interactions, 7 ^> 1. In (|41|1 one can neglect the dependence of the kernel on x and 
y, and g{x) becomes a constant g, which satisfies an equation 



, 80 / ,26 26A , 
^9 = 1 + tan" 1 - + ——- , (45) 



while l|43|l reads 



ttA V A A 2 + 46 2 



-tan- 1 ^. (46) 

7T A 



After some algebra energy is rewritten as 

/ x 71-2 f 4, sin7ra, 12. sin7ra. 9 \ 1 

3 \ 7 7T 7 Z 7T / 7 J 

Using exact solutions, one can analyze demixing instabilitieai2iiii4 for repulsive bose-fermi mixtures. In the absence 
of external potential bose-fermi mixture is stable, if the compressibility matrix 



dfJ.b dfi b 

drib dn / 

djif_ djif 

drib dn / 



(48) 
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FIG. 1: Energy of the ground state is given by E = e(j, a)h 2 N 3 /(2mL 2 ), where 7 = mg/(h 2 n), and a = M/N is the boson 
fraction. When a = 0, system is purely fermionic, and the energy doesn't depend on interactions. When a = 1, the system is 
purely bosonic, and numerically obtained energy coincides with the result oi» 8 . If 7 = 0, bosons and fermions don't interact, 
and e( 7 ,a) = (vr 2 /3)(l - a) 3 . 



is positively defined. Here, rib is the boson density, and rif is the fermion density, fib and [if are the bose and fermi 
chemical potentials, given by 



N 2 ( de de 

!"■ = — ( 3e (7, a ) - 7— + (1 - u)-q^ 

de 



L 2 ^'.-y 1 Ql 

N 2 ( 8e 



(49) 
(50) 



The fact that the matrix (|48|l is positively defined can be checked numerically for any value of a and 7, and proves 
that bose-fermi mixture with the same bose-fermi and bose-bose interactions is stable with respect to demixing for any 
values of bose and fermi densities. We note, that the absence of demixing for one particular value of the density has 
been checked in the original article by Lai and Yangi&. Although an exact solution is available only under conditions 
P)l. small deviations from these should not dramatically change the energy e(7, a). Therefore, we expect the ID 
mixtures to remain stable to demixing in the vicinity of the integrable line J5J for any interaction strength. Recently 
this has been checked numerically in Quantum Monte Carlo studies for a systems of up to 14 atomsi^. 

Note, that prediction oS* about demixing at sufficiently strong interactions in this case is incorrect, since it is based 
on the mean field approximation. Indeed, the demixing condition there reads 



nt < 



gbbfi 2 ^ 2 



For gtf — 9bb and rib = rif it is equivalent to 



mg 
h 2 n 



> 



4.9. 



(51) 



(52) 



Clearly, this condition is incompatible with mean- field approximation, which is valid for 7 < 1. 

For weakly interacting case one can use mean field approximation to calculate energy and chemical potentialsiSiiS 



E = L 



-nl + gn b n f + 



n it n f 



2m 3 



/i 6 = g{n b + n f ), [if = gn b + 



2 ^2 



% 7T 

2m 



■rif. 



(53) 
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For the strong interactions, up to corrections of order I/7 3 

"ftV (n/ + n 6 ) 3 



E 



2m 



A, sin7ra, 12, sin7ra N 
1 - - a + + — (a + , 

7 7T 7 Z 7T 



h it , , 9 ,„ 1 / sin7ra. ,„ 1 

— — (nt + n b y{l + — -16(a H ) + 4a (1 + costto:) + 

2m 07 \ 7r 

1 , sin7ra N / , sin7ra. , , 
+ -j(aH ) 20(aH ) - 8a(l + cos7ra) ), 



2 ^2 



H 7T 



1 



^ 6 = — — (rif + n b ) (1 + — -16(a 
2m - ■ 



sin irct 



1 , 

T 



37 
sin 7ra 



) + 4(a - 1)(1 + cosTra) ) + 



sm 7ra „ , 
20(aH ) + 8(1 - a)(l + cos7ra) ). 



(54) 



(55) 



(56) 



IV. LOCAL DENSITY APPROXIMATION AND COLLECTIVE MODES 

So far our arguments have been limited to the case of periodic boundary conditions without external confinement. 
This is the situation, when the many-body interacting model Q is exactly solvable in the mathematical sense. If 
one adds an external harmonic potential, model is not solvable any more. However, if external potential varies slowly 
enough (precise conditions for the case of bose gas have been formulated in£i), one can safely use local density 
approximation (LDA) to analyze the density profiles and collective modes in a harmonic trap. In the local density 
approximation, one assumes that in slowly varying external harmonic trap chemical potential changes according to 



n , m,uj 2 x 2 



^(0), M°/(z) 



? 1 



(57) 



Let us consider the case when external harmonic confining potential oscillator frequencies are the same for bosons 
and fcrmions. We note, however, that one can also analyze the case when u> b ^ u>f in a similar way. We consider 



LU b =LOf = UJ , 



(58) 



since in this case distribution of the relative boson and fermion densities is controlled only by interactions, and not 
by external potential, since external potential couples only to total density. Eqs. I|57(l for uj b = ojf = ujq imply that 
densities of bosons and fermions in the region where bosons and fermions coexist are governed by 



de 



M ?(0), fi(x) - ${x) = (n f + n b ) 2 — = M °(0) - M ?(0) 



(59) 



One can show, that these equations cannot be simultaneously satisfied for the whole cloud, and the mixture phase 
separates in an external potential given by (|58|l . For both strong and weak interactions bosons and fermions coexist 
in the central part, but the outer sections consist of Fermi gas only. In the weakly interacting limit, this can be 
interpreted as an effect of the Fermi pressure^: while bosons can condense to the center of the trap, Pauli principle 
pushes fermions apart. As interactions get stronger, the relative distribution of bosons and fermions changes, and 
Figs. and 01 contrast the limits of strong and weak interactions. For strong interactions, the fermi density shows 
strong non-monotonous behavior. 

When interactions are small Eqs. (|53(l and l|59(l imply that in the region of coexistence densities are given by 



n° b (x) = n°(0)(l - -), n%x) = n° f (0), for x 2 < x 2 . 

X b 



Outside of the region of coexistence, density of fermions decays as the square root of inverse parabola: 



ng(aO=0, n°(x) 



n° f (0) 



1 - 



X 2 r 1 „ 2 „ 2 

, for x b < x < x 

x f 



(60) 



(61) 



Parameters x / and x b are given by 



2gng(0) 

2 ' 
muiQ 



(W3(0)) 2 
(mw ) 2 



(62) 
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FIG. 2: Densities of bose and fermi gases in weakly interacting regime at zero temperature. Lieb-Liniger parameter in the 
center of a trap is 70 = 0.18, overall number of bosons equals number of fermions. Total density in the center of a trap is taken 
to be 1. 




FIG. 3: Densities of bose and fermi gases in strongly interacting regime at zero temperature. Lieb-Liniger parameter in the 
center of a trap is 70 3> 1, overall number of bosons equals number of fermions. Total density in the center of the trap is taken 
to be 1. 



A typical graph of density distribution for weakly interacting case is shown is shown in Fig. [2 
If effective 70 is much bigger than 1 in the center of a harmonic trap, the total density n°(x) follows Tonks-Girardeau 
density profile 



n°(x)=n°(0)Jl~-j. (63) 




From Eqs. H47|) and (|59|1 distribution of a(x) is controlled by the following equation: 

n°(x) 3 (l + cos (7ra(x))) = n°(0) 3 (l + cos (7ra(0))). (64) 

Since 1 + cos (Tra(x)) is bound and n°(x) goes to near the edges of the cloud, this equation can't be satisfied for all 
x 2 < Xj, which means that only fermions will be present at the edges of the cloud, similarly to weakly interacting 
regime. Density distribution for equal number of bosons and fermions is shown in Fig. [3] The form of the profile is 
universal, as long as 70 3> 1 and the temperature is zero. Evolution of this profile for nonzero temperatures is shown 
in figure ITf! 

Recent experiment*^ demonstrated that collective oscillations of ID gas provide useful information about interac- 
tions in the system. Here we will numerically investigate collective modes of the system, by solving hydrodynamic 
equations of motion. These equations have to be solved with proper boundary conditions at the edge of the bosonic 
and fermionic clouds. Within the region of coexistence of bosons and fermions, such oscillations can be described by 
four hydrodynamic equations^ 

d d 

-7^5n b + —(n b v b ) = 0, (65) 
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d Z Q , T r 1 

m dt + ~dx + + 2 mV ' 



0, 



|fe /+ ^(n /V/ )=0, 
d d 1 



dx 



(66) 
(67) 
(68) 



In certain cases, analytical solutions of hydrodynamic equations are available2^& and provide the frequencies of 
collective modes. When an analytic solution is not available, the "sum rule" approach has been use d2iii2ii2L2& to obtain 
an upper bound for the frequencies of collective excitations. The disadvantage of the latter approach is an ambiguity 
in the choice of multipole operator which excites a particular mode, especially for multicomponent systems^. Here we 
develop an efficient numerical procedure for solving hydrodynamical equations in ID, which doesn't involve additional 
"sum rule" approximation. 

While looking at low amplitude oscillations, it is sufficient to substitute 



n b (x, t) = n°(x) + Sn b (x)e iujt 
rib(x,t)vb(x,t) — n { l{x)8vb{x)e lult 
rif(x, t) = rif{x) + 5nf(x)e 



nf(x,t)vf(x,t) — n°f(x,t)5vf(x)t 



iuit 



fib + V ext b + Ti mv 'b — constl + 5fi b (x)e luJt — constl + {Srib{x)^L + 5rif(x)^^-)e' lujt 
2 drib drib 



1 



Hf + V ex tj + 7. mv f = const2 + 5jjif(x)e 



- const2 + («Wx)§^ + Snt {x)^)e wt 
drib 



dn f ' 



(69) 
(70) 
(71) 
(72) 

(73) 

(74) 
(75) 



Here, ri b (x) and n^(x) are densities obtained within local density approximation. Linearized system of hydrodynamic 
equations can be written as: 



raw 



5n b (x) 
Sn f (x) 



V 



n° b (x) 






n°f(x) 



drib dnj 
drib dnf 



5n b (x) 
Srif(x) 



(76) 



For numerical solutions and boundary conditions it is more convenient to work with independent functions 
S fib(x) , d /j, f (x) . System of equations becomes 



5n b (x) 
6/j,f(x) 



dnb dnf 
djjj_ du f 
dnb dnf 



n° b (x) 






n°Jx) 



6fJ*b(x) 
fix) 



(77) 



Outside of the region of coexistence of bosons and fermions, S/if* satisfies the following equation: 

dnf 



mu'Snf 1 = 



dnf 



■V [n° f ut {x)V8n° f ut ] 



(78) 



All modes can be classified by their parity with respect to x — > — x substitution, and will be investigated by parity- 
dependent numerical procedure. We will consider equations only in the positive half of the cloud. For even modes, 
one may require two additional conditions: 



V8u f (x = 0) = 0, V5u b (x = 0) = 0. 
For odd modes, analogous conditions are 

5n f {x = 0) = 0, 5fj,b(x = 0) = 0. 



(79) 



(80) 



Boundary conditions for fermions at the edge of the bosonic cloud, Xb, correspond to the continuity of Vf and Sfif- 
Continuity of the velocity can be obtained by integrating continuity equation l|67|l in the vicinity of Xb ■ From Eq. 1681) 
it is equivalent to 



VSu° f ut {x = x b + 0) = VSu f (x = x b -0). 



(81) 
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The second condition can be obtained by integrating (|68|l in the vicinity of Xb ■ 

8^° f ut {x = x b + 0) = S^fix = Xfc-0). (82) 

One may see, that these conditions do not imply that 5n°j ut {x = Xb + 0) = 5nf(x = Xf, — 0). This can be easily 
illustrated by the dipolc mode, where 6vf(x) — Svb(x) — const, 5nf = Vriy(x), which is clearly discontinuous for 
profiles shown in Figs. [21 and [3] 

Two additional conditions come from the absence of the bosonic(fermionic) flow at Xb{xf) : 

n°(x)v b (x)\ x - >Xb -o = 0, (83) 
n°(x)vf(x)\ 0. (84) 

Outside of the region of coexistence, the chemical potential and density of fermions are given by /ij u * ~ 

(n°f Ut ) 2 , n°f Ut ~ yl — ( x / x f) 2 i where Xf is the fermionic cloud size. In dimensionless variables u = x/xj, eq. fTHjl can 
be written as 

uj 2 t „ d 2 5u° f ut dSii° f ut 

-^.(l-O-^f— (85) 

For this equation, there exists a general nonzero solution which satisfies (|84|) : 

6tf ut = cos(— arccos — ). (86) 

1 W Xf 

Substituting this into (|81|) - (|82J) . one has to solve eigenmode equations numerically for x < Xb, with five boundary 
conditions (|81|) . (|82(l . l|83|l and (|79|) or (|80J) depending on the parity. These boundary conditions are compatible, only if 
u is an eigenfrequency. Using four of these boundary conditions, the system of two second order differential equations 
can be solved numerically for any uj. To find a numerical solution we choose to leave out condition i|83[) . and check 
later if it is satisfied to identify the eigenfrequencies. 

The most precise way to check (|83f) numerically is based on equations of motion. For even modes, Vb(0) = 0, and 
integrating (|65() from till Xb, one obtains 

Snb(x)dx = (jib(x = Xb)vb(x = Xb) — nb(x = Q)vb(x = 0)1) = 0. (87) 

o «w 

For odd modes, from eq. I|tjtj|) ^&(0) = i\75fib(x = 0)/(mw), and integrating (|65|l from till Xb, one obtains 

r / \ t 1 / / , / . , „. , nb(x = 0)VSiib(x = 0) , . 

dnb(x)ax = (rib{x = XbWblx = Xb) — mix = 0)Vb(x = 0)11 = t. . (88) 

o toj mui 2 

When a numerical solution for 5fib(x),6/j,f(x) is available, conditions H87JI or (|88(l can be checked numerically using 

5n b (x) = — -~ — ^ t, — 5 . (89) 

drij drib drib dnj 

First we apply this numerical procedure for weakly-interacting regime, and the frequencies of collective modes are 
shown in fig. 0] When 70 — > 0, bose and fermi clouds do not interact, and collective modes coincide with purely 
bosonic or fermionic modes, with frequencies^ ojf — nuo^ and tu b — u>q\J n(n + 1) /2. Modes which correspond to 
lu/ujo = l,v3, 2,v6 are shown in fig. 01 As interactions get stronger, bose and fermi clouds get coupled, and all 
the modes except for Kohn dipole mode change their frequency. For Kohn dipole mode, bose and fermi density 
fluctuations are given by Snf = Vn°(a;),(5nt = Vn"(x). In Figs 15171 we show density fluctuations for three other 
modes in the region of coexistence for a particular choice of parameters 70 = 0.394, Xb/xj = 0.6 and equal total 
number of bosons and fermions. Modes for which the frequency goes down due to coupling between bose and fermi 
clouds correspond to the collective excitations with opposite signs in density fluctuations of bose and fermi clouds. In 
TG regime these modes continuously transform into "out of phase" low- lying modes which do not change the total 
density. At weak interactions lowest mode is an "out of phase" dipole excitation, after that comes "in phase" Kohn 
dipole mode (center of mass oscillation), "out of phase" even mode, "in phase" even mode, second "out of phase" odd 
mode. 
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0.2 0.4 0.6 0.8 1 1.2 

FIG. 4: Frequencies of collective excitations in mean field regime versus Lieb-Liniger parameter in the center of a trap 70. 
Total number of bosons equals the number of fermions. Even in mean field regime frequency of "out of phase" oscillations gets 
smaller as interactions get stronger. 
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FIG. 5: Fermi and bose density fluctuations of "out of phase" dipole mode (lj/ujo) 2 = 0.757 for 70 = 0.394, Xb/xf = 0.6. Total 
number of bosons equals the number of fermions. Outside of the region of coexistence of bose and fermi clouds, Snf(x) ~ 
cos(— arccos — ) 



y/l-(.x/x f ) 3 V "0 



Let's consider Tonks-Girardeau regime, when energy is well approximated by (|4T|) . Since dependence of the energy 
on relative boson fraction a(x) is I/7 times smaller than dependence on the total density, the energetic penalty for 
changing relative density of bosons and fermions is small. Thus there should be low-lying modes, which correspond 
to an oscillation of the relative density between bosons and fermions, while total density is kept fixed up to I/7 
corrections. In addition to these low-lying "out of phase" oscillations of bose and fermi clouds, there will be "in 
phase" density modes, which correspond to oscillations of the total density. Since up to 1/7 corrections dependence 
of the energy on total density in TG regime is the same as for free noninteracting fermions, energy of these excitations 
is given by2& to = nu>o, up to small corrections of the order of 1/7. 

When 7—5-00, relative compressibility goes to zero as I/7, so from eq. (|TT|> energy of low- lying modes goes to zero 
as 1/^/to, where 70 is a Lieb-Liniger parameter in the center of a trap. Performing a numerical procedure outlined 
above, one can obtain the dependence of the frequencies of low- lying "out of phase" modes on relative density of 
bosons and fermions. Results of these calculations are shown in fig. [H] and are parameterized by the overall boson 
fraction and 70. It turns out that the lowest lying mode is odd, and after that the parity of collective excitations 
alternates signs. For " out of phase" modes signs of density fluctuations and velocities of boson and fermion clouds 
are opposite. One can easily understand, why does the energy grow, as the boson fraction is decreased: the size of 
the bose cloud shrinks, and the "wavevector" of the corresponding excitation increases, leading to an increase of the 
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arb. units 




FIG. 6: Fermi and bose density fluctuations of "out of phase" breathing mode (w/wo) 2 = 2.51 for 70 = 0.394, Xb/xf = 0.6. 
Total number of bosons equals the number of fermions. Outside of the region of coexistence of bose and fermi clouds, 8rif(x) ~ 




FIG. 7: Fermi and bose density fluctuations of "in phase" breathing mode (lj/luo) 2 = 3.585 for 70 = 0.394, x b /xf — 0.6. Total 
number of bosons equals the number of fermions. Outside of the region of coexistence of bose and fermi clouds, Snf(x) ~ 




FIG. 8: Dependence of the frequency of lowest lying "out of phase" modes for 70 > 1 on overall boson fraction a, where 70 is 
the Lieb-Lininger parameter in the center of a trap. Characteristic scale of "out of phase" oscillations in strongly interacting 
regime is cjo/^/To <S <^o- Total density "in phase" modes tu n = ncoo have much higher frequency for 70 3> 1 and are not shown 
here. 
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frequency. One should note that for very small overall boson fraction 70 3> 1 is not enough to separate energy scales 
for "out of phase" and "in phase" oscillations, and also conditions for applicability of LDA become more stringent. 



V. ZERO-TEMPERATURE CORRELATION FUNCTIONS IN TONKS-GIRARDEAU REGIME 

Calculation of the collective modes in the previous section relies only on the dependence of the energy e(j,a) 
on the densities of bosons and fermions. Collective modes can be used in experiments^Si to check to some extent 
quantitatively the equation of the state of the system^. However, only some part of the information about the 
ground state properties is encoded in the energy: indeed, the energy and collective modes of the strongly interacting 
Lieb-Liniger gas are the same as for the free fermionsSS^, while the correlation functions are dramatically different^. 
Single particle correlation functions can be measured experimentally using Bragg spectroscopy* or time of flight 
measurements^. Generally, it is much harder to calculate the correlation functions compared to the energy from 
Bethe ansatz solution. Most of the progress in this direction has been achieved for the case of strong interactions^. 
Recently, there have been some reports^, where pseudofermionization method has been used to calculate correlation 
functions for spin-^ fermion Hubbard model for the intermediate interaction strengths. In this section, we will analyze 
the correlation functions in the regime of strong interactions, using the factorization of orbital and "spin" degrees of 
freedom similarly to the case of spin-i fermions^i^. Our calculations in this section are performed for the periodic 
boundary conditions, when the many body problem is strictly solvable in the mathematical sense. We will obtain 
a representation of correlation functions through the determinants of some matrices, with the size of these matrices 
scaling linearly with the number of the particles. These determinants can be easily evaluated numerically, and 
provide a straightforward way to study correlation functions quantitatively at all distance scales. This determinant 
representation can be generalized to nonzero temperatures, and results of this generalization will be presented in the 
next section. 



A. Factorization of "spin" and orbital degrees of freedom 

The regime of strong interactions can be investigated in by neglecting fcj compared to A Q , c in (|28 [l -(|29 p . Simplified 
system for spectral parameters is 



-A a + zc/2 x A 
-A Q - ic/2 



= 1, a = {!,..., M}, (90) 



e^=]J -^ + i f j = {l,...,N}. (91) 
- 1 1 —A/3 — ic 2 

(3=1 ' 1 

We see that "spin" part is decoupled from orbital degrees in the Bethe equations. Equation (|90|l for ground state 
"spin" rapidities can be resolved as 



-Kg + ic/2 _ t27rKa / N 



-A Q - ic/2 



a< = {l,...,M}, (92) 



where K a is a set of integer "spin" wave vectors. Since the details of calculations depend on the parity of M and A, 
from now on we will assume that A is even, and M is odd. Ground state corresponds to A a occupying "Fermi sea" 
(—A, A), so from (|92|) ground state "spin" wave vectors are 

n t = {— (M - l)/2 + A/2, A/2, (M - l)/2 + A/2}. (93) 

This choice of "spin" wave vectors will be justified later, in section fvTI From equation Q91|l it follows that ground 
state orbital wave vectors are 

h = {— 7r(A - 1)/L, -tt/L, n/L, n(N - 1)/L}. (94) 
Eq. l(TB|l for F(A, y) simplifies to 
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and "spin" wavefunction l|27|) can be represented as a Slater determinant of M single particle plane waves in "spin" 
space: 

ip = e - J T?I>" det[e*'^ Ki ^]. (96) 

Orbital part of the wavefunction also simplifies into a Slater determinant, since all Yang matrices in (J2J) are equal 
to -1. 

Ground state is written as a product of two Slater determinants, describing orbital and "spin" degrees of freedom: 

$>( Xl ,...,x N ) ~ det[e*^]det[e^ Kl ^]. (97) 

Here x\, % are coordinates of bosons, xm+i, xn are coordinates of fermions, and yi is the order in which the 
particle Xi appears, if the set X \, ...,xn is ordered. In other words, if 

< x Ql < x Q2 < .... < x Qn < L, (98) 
then {y 1 ,...,y N } = {Q- 1 (l),...,Q- 1 (N)}. (99) 

First determinant depends on positions of both bosons and fermions, while the second determinant depends only on 
relative positions of bosons y%, ...,yu- Normalization prefactor will be determined later to give a correct value of the 
density. One can confirm that symmetry properties of wavefunction are as required: transposition of two fermions 
affects only first determinant, therefore wavefunction acquires —1 sign. Transposition of two bosons changes signs of 
both first and second determinants, so wavefunction doesn't change. 

Similar factorization of wavefunction into spin and orbital degrees of freedom has been observed ir4£ for one 
dimensional spin-i Hubbard model. In that case, spin wavefunction is a ground state of spin-i antiferromagnetic 
Heisenberg model, and is much more complicated compared to (|96|l . 

It might seem that "spin" degrees are now independent of orbital degrees, but this is not true, since it is the relative 
position of orbital degrees which determines "spin" coordinates. If one wants to calculate, say, bose-bose correlation 
function, one has to fix position of x\ and x\ and integrate \&(xi,X2, .., xn)^ (x^, X2-~, xjv) over X2, ...,xjv- However, 
there are C^f inequivalent spin distributions, and integration in each subspace l|98|l has to be performed separately. 
For spin-i fermions on a lattice in^ this integration becomes a summation, and it has been done numerically for up to 
32 cites. This summation requires computational resources which scale as an exponential of the number of particles. 
Here we will report a method to perform integrations for a polynomial time, which will allow to go for larger system 
sizes (easily up to 100 on a desktop PC) and study correlation functions much more accurately. 



B. Bose-Bose correlation function 



Let's describe a procedure to calculate bose-bose correlation functions of the model. First, we will use translational 
symmetry of the model to fix the positions of the first particle at points x\ = 0, x' x — £. Instead of writing wavefunction 
as a function of positions of M bosons and N — M fermions, let's introduce a set of N ordered variables 

Z = {0 < zi < z 2 < ... < z N < L}, (100) 

which describe positions of the atoms, without specification of bosonic or fermionic nature of the particle. If any two 
particles exchange their positions, they are described by the same set (|100|l . In addition to (|100|l one has to introduce 
a permutation y which specifies positions of bosons: y\,...,yM are boson positions, and yM+ii •••> UN are fermion 
positions in an auxiliary lattice: z yi = X{. In this new parameterization normalized wavefunction is (normalization will 
be derived later in this subsection) 

y( Zl ,z 2 ,..,z N ;y) = , 1 det[e tfc - Zj ]det[e^ K -^](-lf . (101) 

V ' y/{N -M)\M\L N N M y ' 



Here and later we denote a sign factor 



(-1)"= [] Signivi-yj). (102) 

N>i>j>l 



One should note, that second determinant has a size M x M, and depends only on yi, ...,yM- Dependence of wave 
function on yM+i, JA/v comes only through sign prefactor. For each particular set of y\, yu there are (N — M)\ 



16 



different configurations of j/m+i, UN, for which wave function onfy changes its sign depending on relative positions 
of J/m+i, >~,Vn- 

To calculate correlation function we should be able to calculate a product of wavefunctions at the points 

x\ = 0, x\ — £, x' 2 — x 2 , x' N — x N . (103) 

Let Z be is an ordered set for Xi variables: 

Z = {zi = < z 2 < ... < z N < L}. (104) 

If we denote an ordered set for x\ variables as Z* ', then using l|103|l one can conclude that Z' is obtained from Z by 
removing z\ = 0, inserting an extra coordinate z' d = £, and shifting variables which are to the left of it: 

Z' = {0 < z[ = Z2 < 4 = z 3 < ... < 4_i =z d <z' d = £< z' d+1 = z d+1 < ... <z' N = z N < L}. (105) 

" Spin" states y and y' are connected by 

2/1 = l,2/i = d > 

y'i = y% - 1) for K Vi < d , 

2/i = Vi, for d < 2/i- (106) 

Correlation function can be written as 

p b (0,S,) = M J $(0,x 2 ,..,x N )&(t,X2...,x N )dx 2 ...dx N = 

( (JV - unlt-ll^NM det t e ^ det[e^"^] det[e~^'i] det[e-'»"'^ (107) 

where integration over dzi and summation over y are done subject to constraints (|105fl - H106f) . One can observe now, 
that limits of integration in l|105f) depend only on £ and d. These limits are independent of y, and function under 
integral factorizes into z— dependent and y— dependent parts. Similarly, summation over y doesn't depend on precise 
values of £ or Zi, but the dependence comes through d. Therefore, density matrix can be written as 

1 N 

Ao, o = {N _ mM _ mNm £ m os\d), (los) 

where I(d, £) is a an integral 

l(d,0= [ det[e ikiZ i]det{e- ikiZ 'i]dz 2 ...dz N (109) 



subject to constraints l|105|) . and S b (d) is an expectation value of a translation operator over a symmetrized Slater 
determinant wavefunction: 

S b (d) = ^det[e i ^ K «w]det[e-^%-lf (-If'. (110) 

v 

Normalization can be determined using the following argument: if £ = 0, then only contribution from d — 1 does not 
vanish. One can calculate 1(1, 0) = JVL^ -1 and S b (l) = (N — M)l(M)\N 1 , since these follow from normalizations 
of orbital and "spin" wavefunctions. Since we want p b (0, 0) = M/L, we can fix the normalization prefactor in (|101fl . 

1. Calculation of a many-body integral 

Let's describe the calculation of an integral I(d,^). From now on we will assume that L = 1. First, since ki are 
equidistant wave vectors (|94() . one can use Vandermonde formula to simplify the determinants: 

det[e lfci ^] = e -^(^-i)(^+-+^) det [ e ^(;-i)^] = e -^(N-i)(z 1+ ...+z N ) J"J ^ m _ gfixi,^ j = {i,...,^}, 

Jl<]2 

det[e- iklZ 'i] = e ™(N-i)(z[+-+z' N ) det [ e -i!hr(i-i)*;] = e iK(N-i)(z{+...+z' N ) J I ( e -«hr»J a _ e^* 7 ^), I = {1, N}(111) 
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Using this representation, the fact that z\ = and l|105fl . one can rewrite these N x N determinants as a product of 
[N — 1) x (N — 1) determinant and a prefactor: 



N-l 



det[e lfel ^] = e -*T(^-i)(*i+-+**-0 det[e l27r(; - 1)t! ] ]J ~ 1), i = {1, AT - 1} 

»=i 

iV-l 

det[e lfei ^] = (_i)rf-i e "(^-i)«+*i+-+*«-i) det[e-' i27r(i - 1)t! ] ( e -^ u - e" 12 ^), I = {1, ...,N - 1}, (112) 

z=l 

where we introduced iV — 1 variables of integration tj , so that 

U = z i+1 . (113) 

Factor (— l) d_1 arises since z' d — £, and to write H112|) we changed signs of d — 1 terms in <|llll) . Integration subspace 
is defined as 

{0 < h < ... < t d -! < £ < t d < ... < tpj—i < 1} (114) 

One can extend this subspace as follows: 

T = {0<h,...,t d - 1 <t<t d ,...,t N -i < 1}. (115) 

Indeed, expression under integral doesn't change, when ti < £ and tj < £ change their positions (similarly for U > £ 
and tj > £ ), so this extension just adds prefactor l/((d — 1)1 (N — d)\). Finally, we have 

f_i \d-l p m(N-l)£ f N ~} 
J ( d ^) = u „u AT M / dette^^-^^jdetle-^- 1 )*'] TT (e* 2 ^ - l)( e -^ - e^)^...^!. (116) 
(a — l).(l\ — a). J ti c_T i=1 

At this point we use a trick fromSi, where Toeplitz determinant representation for strongly interacting bose gas was 
derived. Lets expand determinants under integrals using permutation formula for determinants: 

N-l 

j-j ^m-D^ph^u^u _ 1)(e -i2 Wi4 _ e -i^ dtl _ dtN1 ( 118 ) 



i=l 



From summation over P, P' we can go to summation over P, Q, where P' = QP. Also, one can remove constraints 
(|115fl by introducing two functions 

f l (£,t) = (e j27r * - l)(e- i27rt - e- a ^) fort < £, otherwise, 

/ 2 (£,t) = (e l27rt - l)(e~ l2xt - e- l27r «) fort > £, otherwise. (119) 

I(d,£) becomes 

f_-\\d-l„iir(N-l)£ d -! r l r l 

£ £ (-i) Q (Il/ ^-^/^^xn / ^-^)uf^ ti)dtiUm 

( a i M jv aj. PcSiv i QcSjv i l=1 Jo i=d Jo 

If /^jt) and / 2 (£, t) were the same, as in**, expression being summed wouldn't depend on P, and summation over Q 
would give a determinant, with the same elements along diagonals (Toeplitz determinant). In our case, for each given 
P the expression is P— dependent, and the result doesn't have the Toeplitz form. However, introducing additional 
"phase" variable, one can recast the expression as an integral of some Toeplitz determinant. Desired expression has 
the form: 

fl-K . N ~l rl 

I(d,0 = (-l) d - 1 e i *( N -W ^-^%(^(-l)«J] / ^-^(e^f^c^+f^^da), (121) 

Jo 71 q i=1 Jo 
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where Cj is a dummy variable of integration. Integration over ip is analogous to projection of BCS to a state with a 
fixed number of particles. After integration over tp nonzero terms appear, if in the expansion of the product of brackets 
for some d — 1 brackets f 1 is chosen instead of f 2 . If this choice is made at brackets with numbers Pi, Pd-i then 
contribution from such a choice exactly corresponds to a term in (|12L)fl . However, each choice of brackets corresponds 
to (N — d)\(d — 1)! different permutations, and this cancels the same combinatoric factor in the denominator of (|120f) . 
Summation over Q is nothing but a determinant, and finally we have 



I(d,£) = (_l)*-l e *r(W-lK *£. e -Hd-l) V det 

Jo 27r 



co(<p) 



co(y) 



where 



CN-2(V>) 

co(f) 



(122) 



Cj(<p)= / e i ^(e^f 1 (^x) + f(^x))dx 



(123) 



Expression in l|122|) without an integral over tp is a generating function of I(d,£) with the weights e^ v ~ 7r ^ d_1 ' ) , and 
integration over y> extracts a particular term out of this generating function. 

What we achieved in this section is to represent a complicated N — 1 fold integral as an integral over one phase 
variable, which can be done numerically in a polynomial time over N. 



2. Calculation of S b '(d) 

Calculation of S b (d) is very similar in spirit to calculation of the previous subsection. Integration over Xi corresponds 
to summation over j/j, and £ corresponds to d. Final result is a determinant of some matrix. Due to the shift operator 
(|1U6[) this determinant does not have a Toeplitz form, but it is not important for a numerical evaluation. 

We need to calculate 

S b (d) = ^2det[e^ K '^}det[e- l ^ K ^}(^l) y (-l) y ' , (124) 
y 

where Ki is a set (|93|l . Definition of y' according to (|10(jfl can be rewritten as 

Vi = 1, 2/i = d, 

y' i =y i + Si9n{yi - d) -\ i = {2,...,N}, (125) 



where 



Sign prefactor in 11241 can be rewritten as 



Sign(x) = 1, x > 0, 
Sign(x) = -1, x < 0. (126) 



(-i)»(-i)»' - n - n si 9 n (y'i - v'j) = n - d ^= i- 1 ^ 1 - 

i>j i>j i=2 

We see, that (|124ll depends only on yi, ...,um, so from now on we will consider a summation in yi, ...,?/m variables. 
Summation over j/m+i, ■■•! Vn gives a trivial combinatorial prefactor (N — M)l. Furthermore, we can extend possible 
values of j/i, um to yi — Vj,i ^ j, since for such configurations first determinant in l|124|) is 0, and they don't change 
the value of S b (d) : 

y = {i<y2,ys,-,VM<N}. (128) 

Lets use the fact that k, is a set of equidistant numbers (|93|l , and rewrite determinants using Vandermonde formula, 
similar to (|lllfl : 

det[e l ^ KiW ] = e ^(-(^- 1 )/2+JV/2)(i+...+aA/) det[e l ^ ( '~ 1)w ] = 
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e ^(~(M-l)/2+N/2)(l + ...+y M ) Yl ( e i2 H-V32 - e^^ 1 ), / = {1, M}, 

j'Kj'2 

det[e"^ Kl ^] = e~ 1 ^ (-(A-f-i)/2+JV/2)(d+...+ y ; f ) d e t[e^ = 

e -i%-(-(M-l)/2+N/2){d+...+y' M ) J-J (gt^w^a _ gi^y'n)^ _ {1,...^}. 



(129) 



For simplicity of notations later, lets introduce £j = y»+x , t£ = t/.- + 1 ,i = {l,...,M— 1}. Analogously to (|112fl , we extract 
a determinant of (M — 1) x (M — 1) matrix out of Vandermonde product: 



M-l 



det[e^ K ^] = e ^(-(M-i)/2+JV/2)(i+t 1 ...+tM- 1 ) det^^-i)*;] [] (gifr* _ e** 1 ), i = {1, M - 1}, 



1=1 

M-l 



det[e" l * K ^] = e -m-(M-i)/2+N/2)(d + t' 1 ...+t' M _ 1 ) det [ e -i^(l-i)tJ] JJ ( e -<frt, _ e -^<*), / = {i, - 1}. (130) 



i=l 



At this point we need to represent the subspace of summation (|128[) as a sum over M inequivalent partitions, similar 
to representation (|108fl : 



A I 



b ( r _i)!(M-r)! ( ' j ' 



where S b (d,r) is a result of summation in the T r subspace: 



T r = {1 < tl, ...,*r-l < d < t r , ...,t M -l < N}. 



(131) 



(132) 



Note, that S b (d,r) = for r > d, since in this case two of t\, ...,t r -% should coincide, and wavefunction becomes 0. 
Calculation of S b (d, r) is very similar to calculation of /(£, d). Let's expand the determinants H13Ufl using permutations: 

S\d,r) = {-l) d - x e i2 S'^ M - x y 2+N ' 2 ^-^ Y, E (-l) P (-l) P 'x 



PCSm-i P'GSm-i 



-—1 d 



TT( V e^^-^-^-^-^ie*^** - e i ^)(e-^ ( * I - 1) - e" 4 ^)) x 



i=l ti=X 



M-l N 

Il( E e''^ ((p * -1) - (P4 - 1))t4 (e^** -e^)(e _i ^' -< '^' ? 

i=r t i= d+l 



)) 



(133) 



From summation over P, P' we can go to summation over P, Q, where P' = QP. Also, one can analytically perform 
summation over ti in each of the brackets, since it is a combination of geometrical progressions(this is analogous to 
integration over ti variables in previous subsection): 

r-1 M-l 

S b {djr) = { _lY-l e m-{U-l)/2 + N/2) { r-d) £ £ (-1)«JJ C^Q^Pi) J] C 2 (d,g P „P), (134) 

PGSm-i QcSm-i 



i=l 



where 



t=d 



t=7V 



2 (d,j,l)= e l ^W- i )*(e i ^ d -e I ^)(e- l ^*-e- l ^ d ) 



(135) 



t=d+l 



are independent of r. At this point, we can use the "phase" variable integration trick to get rid of summation over P, 
and then represent summation over Q as a determinant: 



^_ e -i(r-DV det 
2vr 



5 b (d,r) = (r- l)!(M-r)!(-l) d - 1 e i ^(-( M - 1 )/ 2+JV / 2 )^ x 

c(V>,l,l) c(^,2,l) ... c(V,M- 1,1) 
c(^l,2) c(V>,2,2) ... C (V,M-1,2) 



c(^,l,M-l) c(^,2,M-l) ... c(V>,Af-l,M-l) 



(136) 
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where 

c(ij,jj) = e^c 1 (d,j,l)+c 2 (dj,l). (137) 
We can analytically perform summation over r in (|131f) . since the determinant and c(ip,j,l) are independent of r : 

S \d) = J2 tSS^S'Wr) = (N-M)\(M-l)\(-l) d - 1 e- i ^ M -»/ 2 + N ^ d -^ x 



^ (r-l)\(M-r)\ 

^ ^(V e i(^(-(M-l)/2+JV/2)-V)(r-l)) det 



j , M 



c(V,l,l) c(V,2,l) 
c(V>,l,2) cfoM,2) 



c(V,M- 1,1) 
c(^,M- 1,2) 



(138) 



c(V>,l,M-l) e(^,2,M-l) ... c(^,M-l,M-l) 
Expansion of the determinant (|138|l in a series over e 1 ^ has terms up to e l< - M_1 ^ : 

A/-1 

det(^) = Y, 

n=0 

Summation over r and integration over t/> lead to 

2tt t / M ,,27r . , Af-1 M 

^(^ e '(f(-(M-l)/2+JV/2)-*)(r-l)) det ^) = / 2HL Y ^/ ne i(Tr(-(JW-l)/2+JV/2)-V)(r-l)+^n = 



(139) 



2tt 



A/-1 



n— r=l 

Y /ne^ ( " (M ' 1)/2+7V/2)n = det(^(-(Af-l)/2 + 7V/2)). (140) 

n=0 

Finally, if we introduce a notation tp = 2vr(-(A/ - l)/2 + N/2)/N, 

S b (d) = (N — M)!(M - i)!(_i)rf-i e -^(-(M-i)/2+^/2)(d-i) x 



det 



c(Vto,M) c(^,2,l) 
c(V>o,l,2) c(V-o,2,2) 



c(V>o,M - 1,1) 
c(^),M- 1,2) 



c(^o,l,M-l) c(^o,2,M-l) ... c(V>o,M-l,M-l) 



(141) 



C. Fermi- Fermi correlation function 



Calculation of fermionic correlation function closely reminds the calculation of Bose-Bose correlation function, so 
we will be sufficiently sketchy in our derivation. First, one splits integration into integration over orbital coordinates 
Zj from the set 



Z = {0 < z x < z 2 < ... < z N < L}, 



(142) 



and summation over "spin" variables. Integration over orbital variables is absolutely identical to the Bose-Bose case, 
the difference comes only from "spin" part S'(d) : 



p f (0,0 = (N-M) / V(xi 1 X2,...,Q)&(xi,X2...,Z)da; 1 ...dx N - 1 



N 

EE 

d=l y 



(-1)»(-1)* 



(N-M- \)\M\L N N M 



det[e*^]det[e l * Kl »]det[e- lfcl ^]det[e- l * K '^] dz 2 ...dz N = 



N 



(N-M- 1)\M\L N N M 

where I(d,£) is given by Ijl22|l . and 

Sf(d) = J]det[e l ^' t *^]det[e"^ Kl ^](-l) y (-l) ? ''. 



J2Hd,0S f (d), (143) 



(144) 
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In (|144fl y' and y are related by 

vn = l,y' N = d, 



, Sign( yi - d) - 1 . 

y l =y l + 2 ,t = {l,..., N-l}. (145) 



Similarly to 1127|l sign prefactor can be rewritten as 



N-l 



(-l)y(-l)y' =l[Sign(y l -y J )l[Sign(y>-y' 3 ) = (-l) N - 1 ]J Sign(d - Vj ) = (-l)^ 1 . (146) 

i>j i>j j=l 

We see , that (114411 depends only on yi, j/Mj so from now on we will consider a summation in yi, tjm variables. 
Summation over J/m+Ij ■■■iVn-i gives a trivial combinatorial prefactor (N — M — 1)!. Furthermore, we can extend 
possible values of yi, ...,J/m to j/j = 7^ j, since for such configurations first determinant in 1144(1 is 0, and they 
don't change the value of S-f(d) : 

y = {2<yi,y 2 ,...,yM<N}. (147) 

We can to represent the subspace of summation (I147fl as a sum of M + 1 inequivalent partitions, similar to 
representation l|131l) : 

f ™ (JV-M-1)!(M)! . 
5 W = lj (r-l)!(M-r+l)! g (< *' r) ' (W8) 

where S* (d, r) is a result of summation in the T r subspace: 

T r = {2 < ti, ..,t r _l < d < t r , ...,t M < N}. (149) 
Product of two determinants in 11441) is rewritten as 

det[e^ Kl W]det[e- l * K '^] = jW-W-V/t+WK'- 1 ) det[e'^ det^^H'- 1 ^], I = {l,...,M}. (150) 
We can expand the determinants (|15()|) using permutations: 

Sf(d,r) = (-l) d - 1 e i ^-( M - 1 V 2+N W r -V E (-l) P (-l) P 'x 

PCSm P'CSm 
r-1 d M N 

JJ(£j e i^(W-i)*«-(fl-i)(**-i)))JJ( £ e^CC^-D-^-i))**). (151) 

i=l t i= 2 i=r tj=d+l 

From summation over P, P' we can go to summation over P, Q, where P' = QP. Also , one can analytically perform 
summation over ti in each of the brackets, since it is a geometrical progression. 

r-1 M 

S f(d,r) = (-ir-^-( M -W +N / 2 ^ £ E (-1) Q 114(^^-^)114^^0^ (152) 

PCSm QcSm »=1 i=r 

where 



c}(d,j,0 = e^( I - 1 )X; 



t=2 
t—N 



(d,j,l)= E e^Ci-O* (153) 



4 

t=<2+l 



are independent of r. At this point, we can use the "phase" variable integration trick to get rid of summation over P, 
and then represent summation over Q as a determinant: 



Sf(d,r) = ( r -l)!(M-r + l)!(-l) d -V^(-( M - 1 )/ 2 +™'-- 1 ) x 

c/^,1,1) c/^,2,1) ... C/ (^,M,1) 
c/^,1,2) c/^,2,2) ... C/ (V,M,2) 



27r # e -i(r-l)V det 
2tt 



C/ ftM,M) C/ (^,2,M) ... c f (ip,M,M) 



(154) 
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where 



c f (^j,l)=e^c 1 f (d,j,l)+c}(d,j,l). 



(155) 



We can analytically perform summation over r in l|148|l . since the form of the determinant and Cf(tp,j,l) are inde- 
pendent of r, and r— dependent combinatorial prefactor cancels: 



M+l 

S (r-l)!(M 

27r e 4( T? ("(M-l)/2+iV/2)-0)(M+l) _ 1 



2lT e i(^(-(M-l)/2+JV/2)-V) _ i 



det 



1) 



C/ (^,M,1) 
c/(^M,2) 



C/ (^,1,M) C/ (^,2,M) ... c f (^,M,M) 



(156) 



Analogously to the case of bosons, integration over -0 is equivalent to substitution ipo = 2tt(— (M — l)/2 + N/2)/N to 
the determinant, and the final expression is 



Sf(d) = (N — M - l)!(M)!(-l) d - 1 det 



c/#o,M) c/(V»oj2,1) 
c/^0,1,2) 0/(^0,2,2) 



c/(^,M,l) 
c f (ij ,M,2) 



_c f {i) Q ,l,M) c f (ipo,2,M) ... c/(^,M,M) 



(157) 



D. Numerical evaluation of correlation functions and Luttinger parameters 

Using results of the previous sections, one can evaluate correlation functions on a ring numerically and extract both 
long-range and short range behavior of correlation functions. Calculation of all determinants requires polynomial 
time in their size, and systems of up to N = 100 atoms can be easily investigated on a desktop PC. Fourier transform 
of correlation function is an occupation number n(k), which can be measured directly in time-of- flight experiments 8 
or using Bragg spectroscopj^. Recently, long distance correlation functions of the model under consideration have 
been investigated based on conformal field theory (CFT) arguments 1 "^. Our determinant representations for strongly 
interacting mixture can be used to obtain these correlation functions at all distances, and compare their large distance 
asymptotic behavior with predictions of CFT. 

In figlHlwe show numerically evaluated Bose-Bose correlation function for M = 15, N — 30. Since we used periodic 
boundary conditions, correlation function is periodic in £. To extract universal long-distance correlation functions 
from our calculation, one has to fit the numerical results using general Luttinger liquid asymptotic behavior. In the 
thermodynamic limit long range behavior is 

p b (0,£)Hr 1/(2i H (158) 

where Kb is a bosonic Luttinger Liquid parameter. This formula is valid, if £ is bigger then any non-universal short- 
range scale of the model. In our case, such short-range scale is given by the interbosonic distance, which is L/M. For 
a finite size system, general arguments of conformal invariance^*^ imply that correlation function has the form 

P b (0,O- , . J 1/(2K y (159) 

We fitted numerically obtained correlation functions with 1159(1 . and results coincide with the formula 

(a — l) z — 1 

obtained ini£ based on CFT arguments. One can see subleading oscillations in the numerical evaluation, but their 
quantitative analysis would require more numerical effort. Fourier transform of p b (0,£,) is a monotonously decreasing 
function, which has a singularity at k = 0, governed by Luttinger liquid parameter Kb : 

n\k) ~ | fc| -1+1/(2*6) for k o ( 161 ) 

Fermionic correlation functions can also be obtained using the results of the previous section, and space dependence 
of a typical correlation function is presented in figure 1101 Oscillations are reminiscent of Friedel oscillations of the 
ideal fermi gas. Their large distance decay is controlled by Luttinger liquid behavior. 



23 



P b (0,§ 




4 ^ ^ ^ ^ ^ \ 
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FIG. 9: Normalized bose-bose correlation function on a circle as a function of the distance £ (here rib is bose density). Due 
to periodic boundary conditions correlation function is periodic with a period L, where L is the size of the system. Numerical 
evaluation is done for M — 15, N — 30. 




FIG. 10: Normalized fermi-fermi correlation function on a circle as a function of the distance £ (here n/ is fermi density and 
L is the size of the system). Oscillations are reminiscent of Friedel oscillations of the ideal fermi gas, and their large distance 
decay is controlled by Luttinger liquid behavior. Numerical evaluation is done for M — 51, N = 100. 

One can investigate Fourier transform of the correlation function, which is an occupation number, and results for 
different boson fractions are shown in figs. II 11131 In figure ^2 densities of bosons and fcrmions are almost equal. Fermi 
step at kf gets smeared out by interactions, but relative change of occupation number as kf is crossed is significant. 
As boson fraction is decreasing, the discontinuity appears at kf + 2k{,, and it gets stronger as M/N decreases (see figs. 
I12I13|1 . The presence of this discontinuity has been predicted inii, based on CFT arguments, and here we quantify 
the strength of the effect. One should note, that discontinuity at kf + 2kb is a direct signature of the interactions and 
its detection can serve as an unambiguous verification of our theory. 

VI. LOW TEMPERATURE BEHAVIOR IN TONKS-GIRARDEAU REGIME 

In the previous sections we considered density profiles and developed an algorithm to calculate the correlation 
functions of the ground state of the bose-fermi hamiltonian in the strongly interacting regime. An important 
question, which is very relevant experimentally, is the effect of finite temperatures. In principle, one can use techniques 
of the thermodynamic Bethe ansatz— to obtain free energy at nonzero temperatures as the function of interaction 
strength and densities. Combined with local density approximation, it can be used to calculate density profiles 
for any interaction strength. In this section we will limit our discussion to effects of small nonzero temperatures 
T <C Ef = (ntm) 2 / (2m) only for strongly interacting regime. We will show the evolution of the density profile (see 
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FIG. 11: Fourier transform of the Fermi-Fermi correlation function for M — 51, N = 100. Fermi step at kf gets smeared out 
by interactions, but relative change of occupation number as kf is crossed is significant. 




FIG. 12: Fourier transform of the Fermi-Fermi correlation function for M — 31, iV = 100. Fermi step at kf gets smeared out 
by interactions, and additional discontinuity appears at kf + 2kb- 
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FIG. 13: Fourier transform of the Fermi-Fermi correlation function for M — 11, N = 100. Discontinuity at kf + 2kb gets 
stronger as M/N decreases. 
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fig. I14|) in a harmonic trap and calculate the correlation functions under periodic boundary conditions. The effect 
of nonzero temperatures on correlation functions is particularly interesting for strongly interacting multicomponcnt 
systems (as has been emphasized for the case of bose-bose and fermi-fermi mixtures ir4£) , due to considerable change 
of the momentum distribution in the very narrow range of the temperatures of the order of Ef /j. For the case of bose- 
bose or fermi-fermi mixture it was possibles^ to obtain correlation functions only in the two limiting cases T <C Ef /j 
and Ef/"/ <C T <C Ef. For bose-fermi mixture, we are able to calculate correlation functions for any ratio between 
Ef /j and T <C Ef (see fig. I15f) . By adding an imaginary part to T, the procedure presented in this section can be 
also easily generalized for non equal time correlations. 



A. Low energy excitations in Tonks-Girardeau regime. 

As has been discussed in section ITTT1 for 7 ^> 1 there are two energy scales in the problem: the first energy scale 
is the fermi energy of orbital motion Ef = (irhn) 2 / (2m) , while the second is the the "spin wave" (relative density 
oscillation) energy Ef /j. The second energy scale is present only in strongly interacting multicomponent systems, as 
has been emphasized carliei4Siil. Density profiles and correlation functions we have considered earlier are valid in the 
regime, when temperature is smaller than both of these energy scales: 

T^LEf/KLEj. (162) 

However, interesting phenomena— 4L48,4.9 can ^ e ana jy zec i j n t ne "spin disordered" regime, when 

£//7<CT<c£/. (163) 

This regime has attracted lots of attention recently in the context of electrons in Id quantum wires^lM^. In "spin 
disordered" regime, "spin" degrees of freedom are completely disordered, while orbital degrees are not affected much. 
From the point of view of orbital degrees, this is still a low-temperature regime, since T <C Ef. The energy of the 
system doesn't change too much, while momentum distribution changes dramatically as temperature changes from to 
the order of several Ef/j. "Spin disordered" regime exists only for multicomponent systems and a crossover from true 
ground state to "spin disordered" regime provides a unique opportunity to study the effects of low temperatures on a 
highly correlated strongly interacting system. "Spin disordered" limit is likely to be reached first in the experiments, 
and a significant change of the density profile and of the momentum distribution as regime (|162|l is reached can be 
used as a way to calibrate the temperatures much smaller than Ef. 

Only two limiting cases (|162f> - l|163fl have been investigated for spin-i fcrmion and boson mixtures, since in these 
cases "spin" wavefunctions are related to eigenstates of spin-i Heisenberg hamiltonian, and have a complicated 
structure. In the case of bose-fermi mixture, "spin" wavefunctions correspond to noninteracting fermionized single- 
spin excitations, and one can calculate correlation functions in the whole low-temperature limit, investigating crossover 
from true ground state to "spin disordered" limit: 



Ef/^T^Ef. (164) 

In the following calculations, we will neglect the influence of nonzero temperature on orbital degrees, and will always 
assume that orbital degrees are not excited. This assumption will affect the results only at distances, at which the 
correlation functions are already very small due to effects of spin excitations. 

In the zeroth order in I/7 expansion, energies of all spin states are degenerate, and solutions of Bethe equations 
are given by 

{1,...,M}, (165) 
{1,...,N}. (166) 

In the next order in 1/7 expansion, both kj and A$ acquire corrections of the order of 1/7. Since energy depends only 
on p(k), we need to calculate corrections to p(k) in the leading order. According to l|35fl . to calculate 1/7 correction 
to p(k), one can use Aj in the zeroth order, given by l|165fl : 

1 M A 

1=1 1 



-Ar 



■ic/2 



-A a - ic/2 

M 



N 



1, a 



Ar 



0-- 



ic/2 

\ -A fi -~ic~/2 



,3 
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is independent of k in the first order of I/7 expansion. If we define "spin" wave vectors according to 



-A Q + ic/2 



i2TTK a /N 



a = {!,... ,M}, (168) 



-A Q - ic/2 

energy of the state with " spin" wave vectors Hi in 1/7 order is given by 

£(7, *) = y ^ - 1 £(1 - cos f=p)). (169) 

' i=i 

Allowed values for "spin" wave vectors are 

K = {n t C {1, N}, m < Kj for i < j.} (170) 

The number of "spin" excitations (we will call them magnons from now on) is fixed to be the number of bosons, 
and different "spin" wave vectors cannot coincide. Hence, magnons have a fermionic statistics. The effect of nonzero 
temperatures is to average the correlations over the different sets of possible m from Q17()(l . 

According to l|169|l in the first order in 1/7 expansion magnons do not interact with each other, and the total energy 
is the sum of separate magnon energies. Magnon energy spectrum is 

. . 4tt 2 N 2 , 2nn . 4E f , 2nn 

= W (C0S — - 1} = -37 (C ° S — - 1} - 

Lowest state corresponds to k = N/2, and as the number of magnons increases, "spin" wave vectors k near N/2 start 
being occupied - l|171fl proves the choice (|93|l for the true ground state at zero temperature. 



B. Density profiles 

In this subsection we will analyze the behavior of the strongly interacting mixture in a harmonic trap at low 
temperatures. Similarly to section llVl we consider the case 

UJf, = LOf = LOq. (172) 

According to (|59|l . within the region of the coexistence densities are governed by equations 

fi {x ) + H^tl = ^(o), ^(x) - fi(x) = ^(0) - ^0(0). (173) 
Similarly to the case of T — 0, total density is given by 1631) : 



n Q (x)=n°(0)^l~^, (174) 

and has a weak temperature dependence. On the other hand, relative density is controlled by solutions of the second 
equation l|173f) . and its dependence on temperature is quite strong. It turns out, that in strongly interacting regime 
Mb ~ M/ can be easily calculated using formulas from the previous subsection, fib — fif is the change of the free 
energy, when one boson is added and one fermion is removed from the mixture. On the language of the magnons this 
corresponds to an addition of one magnon. Therefore, one obtains 

Mb^M/^Mm, (175) 

where fj, m is the chemical potential of the magnons with energy spectrum ()171f) . As has been noted earlier, magnons 
obey fermionic statistics (only one magnon can occupy each state) and do not interact, so one can use Fermi distri- 
bution for their occupation number. Chemical potential for magnons \i m as a function of a and T can be obtained 
numerically from the normalization condition for the total number of magnons, which reads 

f 27r 1 dk , s 

a = / is (176) 

Jo e T(^r( cosfc - 1 )-A' m ) + 1 27r 

After that, one can use LDA to obtain the density profiles. In fig. ^]we show the density of fermions for the case, 
when total number of bosons equals total number of fermions. One sees, that density profile changes considerably at 
the temperatures of the order of E'j/jo, where E® and 70 are the Fermi energy and Lieb-Liniger parameter in the 
center of the trap. For E*j/j <C T <C E® boson fraction a is uniform along the trap. As temperature is lowered, 
more bosons condense towards the center of the trap, and fermionic density behaves non-monotonously as a function 
of the distance form the center of the trap. 
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FIG. 14: Evolution of the fermionic density profile in the strongly interacting regime as the function of temperature. Four 
graphs correspond to temperatures O.OLE//70, 0.2Ef/jQ, £7/70 and lO-E^/70. Here 70 is the Lieb-Liniger parameter in the 
center of the trap, Xf is the size of the fermionic trap, and Et = (-Khrio) 2 / (2m) , where no is the total density in the center 
of the trap. Overall number of bosons equals the number of fermions. Non monotonous behavior of the fermi density profile 
persists up to T ~ Ef/jo. The total density profile doesn't change considerably in this range of the temperatures, and is given 
by n(x) = n° y/1 — x 2 /x 2 . 



C. Fermi- Fermi correlations 



From now on we will consider the periodic boundary conditions, when the many body problem is strictly solvable in 
the mathematical sense. We will first describe the calculation of fermi correlations, since it is simpler than calculation 
of Bose correlations. To calculate temperature averaged correlation functions, we should be able to calculate 

f/n , m\ S K , c ^e~Si e(K * )/T (N- A/) f *(ki,...,k A /;xi,x 2 ,...,0)^^ 
9 s (0, £; T) = „ . (177) 

Denominator in i|177|) is a partition function of noninteracting fermions in a micro canonical ensemble. It can be 
written as 



KiCK 



de N 



Z = £ e-2> (Ki)/T = / e- M9 ^l[(l + e ie e^/ T ) (178) 



K = l 



Numerator can be simplified using the factorization of "spin" and orbital parts, similarly to (|143fl : 

T) = Z (N - M - 1VMW ^ ^ e-^« K >V T Sf(K 1 ,..., KM ;d)I(d,£;K U ..., KM ). (179) 

Here S^(ki, km! d) is an expression (|144fl for an arbitrary choice of Ki belonging to 1)170(1 : 

S f (Ki,...,K M ;d) = ^det[e i ^ K *w]det[e-^ Ki ^](-l) s '(-l) J ''. (180) 
y 

I{d, £; K\, Kj\/) is an integral l|109|) . which dependence on K\, % comes only through boundary conditions (|166f) . 
If K i m °d N = D, where D = {1, TV— 1}, then the set of ki which minimizes kinetic energy is uniquely defined: 

A-N/2 + D/N), ^(-N/2 + 1 + D/N), - 1 + D/N)}. (181) 

Li Li Li 

If D — 0, then there are two degenerate sets of ki, and each of them should be taken with a weight 1/2. Taking this 
into account, I{d, £; Ki, km) can be expressed as 

J(d, C; k 1s km) = J(d, 0^(1- - ^ ^)e" (D - f ^ , (182) 

D=0 i=l 
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where 

S N (x) 

5n(x) can be represented as a Fourier sum, 



1 if x mod N = 0, 
otherwise 



N-l 



5n ^ = jf E e ^ px - 



(183) 



(184) 



p=0 



Taking this into account, correlation function (|179|) is rewritten as 



/(0,£;T) 



Z (N - M - l)!M!L Jv iV M + 1 ^ 

v ; d=l 



N N - , , N-l 

rE^-OE(l - 2^)e-(«-*)V« ^ e^S^T), (185) 



D=0 



p=0 



where 



Sf(d;p-,T)= J2 e-^^ + ^^Sf(K U ..., KM ;d). 



(186) 



KiC-R" 



Calculation of S^(kx, km; d) closely reminds a calculation of S^(d) in section IY CI so we will present only a brief 
derivation. 



of/ n (JV-M-l)l(M)! g/ , 

5 J ( K i,..., KM ;rf) = 2^ ( r -i)!(M-r + D! ( M ' ' 



M+l 

^ '(r-l)!(M-r + l)! 
where S*(ki, km', d, r) is a product of two determinants: 

Sf(K U ... ) K M ;d,r) = (-l) d - 1 J2 E (- 1 ) P (- 1 ) P ' 



(187) 



PCSm P'CSm 



r-l d M N 

i—l t i= 2 i=r t 4 =d+l 



(188) 



From summation over P, P' we can go to summation over P, Q, where P' = QP. Also, one can analytically perform 
summation over ti in each of the brackets, since it is a geometrical progression. 



M 



ki, km 



(189) 



PCSm QCS m 



where 



g}(d,3,l)=e^ l J2e l ^- l) \ 

t=2 

gj(d,j,l)= £ 

t=d+l 



(190) 



are independent of r. We can use the "phase" variable integration trick to get rid of summation over P, and then 
represent summation over Q as a determinant: 



2tt' 



e -i(r-l)+ det 



km; d, r) = (r - 1)!(M - r + l^-l)^ 1 x 

5 / (-0,ki,ki) g / (V', Kl,K 2 ) ... g f (lp,Kl,K M ) 
5 / (-0,K 2 ,Kl) 3 / (V', K 2 , K 2 ) ... g f (^,K 2 ,K M ) 

jf(lp, KM,«l) 9 f (lp,KM,K2) ... g f (lp,KM,K M ) 



(191) 
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where 



gf(4 ) ,j,l)=e i ^g}(d,j,l)+gf(d,j,l). 



(192) 



We can analytically perform summation over r, since the form of the determinant and g^(ip,j, I) are independent of r, 
and combinatorial prefactor cancels in i|189|) . Similarly to l|141|) we represent summation over r and integration over 
ip as a substitution tpo — 0, and obtain the following result: 



S f (ki, ...,K M ;d) = (N- M - l)!Af!(-l) d-1 dct 



3 / (0, «i,/ci) 3 / (0,k 1 ,k 2 ) 

3 / (0, K 2 ,Ki) 3 / (0,K 2 ,K 2 ) 



5^(0, «i, «m) 
g / (0, k 2 ,k m ) 



g f (0,KM,Kx) 5 / (0,k m ,k 2 ) ... g f (0,n Ml K M ) 



(193) 



To calculate S' / (d; p; T) we have to sum ()193l) for different choices of Ki with «j dependent prefactor. One can take 
these prefactors into by multiplying each row in (|193[1 by 



f( Ki ) = e-i^w-p^+^m^ 

since only one term from each row appears in the expansion of the determinant: 

S f {d; P ;T) = (N — M — l)!Af!(— l) d ~ 
/(«i)5 / (0, ki,ki) f(K 1 )g f (0,K 1 ,K 2 ) ... f{Ki)g f (0, «i, km) 

/(«2)ff / (0, K2,«l) f(K 2 )g f (0,K 2 ,K 2 ) ... f(K 2 )g f (0,K 2 ,K M ) 



(194) 



A' 



E ••• E det 



/(km)s^(O, «m,«i) /(^m)s' / (0, k m ,k 2 ) ... /(km)S' / (0, k m , km) 



(195) 



Summations over «j in l|195|l can be performed analytically, since each choice of Kj is a term in the expansion of the 
Fredholm determinant 50 . The desired expression has the form: 



S f (d;p;T) = (N — M — 1)!M!(— l) d_1 



e i(N-M)B x 



det 



e <fl + /(l)^(0,l,l) /(^(O.M) 
e*> + /(2)/(0,2,2 



2tt 

f(i)gf(0,l,N) 
f(2)gf(0,2,N) 



/(%'((), 2,1) 

f(N)gf\o,N,l) f(N)gf(0,N,2) ... e ie + / (N)gf (0,N,N) 



(196) 



Integration over extracts terms from the determinant which have e l { N - M ) e dependence. Such terms appear, when 
N — M e 10 elements in the expansion of the determinant are taken along the diagonal. If e are chosen in the rows 
except for Ki, km, then contribution from such choice of e l8 is a minor which equals /(ki).../(km)5 i ^(^i, d). 
Thus evaluation of the prefactor in the e K N -M)6 dependence of the determinant corresponds to summation of 
f(ni)...f(nM)Sf(Ki, KM',d) over possible sets of K{. 

Finally, substituting Ijl96(l into Ijl85|l . one can evaluate numerically fermi-fermi correlation functions for any tem- 
perature and ratio between boson and fermion density in low temperature limit. 

In fig. [E|we show numerically evaluated fermi-fermi correlation function for M = 15, N = 30 and several temper- 



atures, ranging from T = to T = 5^f . At this low temperature region fermi-fermi correlation function changes 
considerably due to transition from true ground state to "spin disordered" regime. In "spin disordered" regime fermi 
singularity at kf gets completely smeared out by thermal "spin" excitations. 



D. Bose-Bose correlation function 



Bose-Bose correlation functions also change as T goes up. However, since for T = n b (k) doesn't have any 
interesting structure except for singularity at k = 0, the effects of nonzero temperatures will not be as dramatic as for 
fermi correlations. We present here the results mainly for the sake of completeness. Calculations in this subsection 
are similar to what has been done in the previous subsection. Correlation function can be written as 

bfn t _ A E^cif e ^^' e{Ki)/T M J ^(Ki, k m ;0,x 2 , ...,x n )^(ki, ...,K M ;^x 2 ...,x N )dx 2 ...dx N 
P(y,€;l) = — , w „ • (!»') 
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FIG. 15: Fourier transform of the Fermi-Fermi correlation function for M = 15, N = 30. Four graphs correspond to tempera- 
tures 0, 0.1-g- 2 -, 0.5-^-, 5-^-*-. In the range of the temperatures ~ -E//7 <C Ef fermi correlation function changes considerably 
due to transition from true ground state to "spin disordered" regime. In "spin disordered" fermi singularity at fc/ gets completely 
smeared out by thermal "spin" excitations. 



Similarly to 1|185J) . this can be written as 



Z(N- M)\(M - l)\L N N M + x ^ 



N N x , , N-l 



D=0 



p=0 



where 



S\d;p;T) = E e-^^ + « K ^S\ Kl ,..., KM ;d). 



KiCK 



Here S b (ni, km; d) is an expression 1)124(1 for an arbitrary choice of Ki belonging to 1(170|1 : 
S b (K 1 ,...,K M ;d) = ^det[e i #*W]det[e^*^](-l)»(-l)*'. 



(199) 



(200) 



Similarly to l|131fl . it can be written as 



M 



S b (Ki,...,K M ;d) = } 



{N-M)\{M-l)\ obi 
^ ( r _l).(M-r)! S (^-^M;d,r), 

where S b (Ki, km', d, r) is a result of the summation of H200|) in the following subspace: 

{1 < V2, -,y r <d< y r+1 , ...,y M < N}. 
We can expand determinants of l|200(l using permutations: 



(201) 



(202) 



S b ( Kl ,... lKM ;d,r) = (-l) d - 1 E E (-l) P (-l) P 'x 

PCS M P'CSm 
r d M N 

i=2 yt=2 i=r+l yi =d+l 



(203) 
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From summation over P, P' we can go to summation over P, Q, where P' = QP. Also , one can analytically perform 
summation over j/j in each of the brackets, since it is a geometrical progression. Compared to the case of fermions, 
there are 3 types of the brackets: 



M 



S f {Ki,...,K M ;d,r) = (-l) d - 1 (-^ Q 9bid,K Q ^, KPt )]Jgl{d, KQpi , KPi )l[gl(d, KQPi ,K Pi ), (204) 

PCS M QCSm 



where 



i=1 



g° b (d,j,l) = e^ 1 - ld \ 
gl(d,j,l)=e i % l J2e i ^-W j 

t=2 

](dj,l) = £ e**W-0* 

t=d+l 



9b ( 



(205) 



We can use "phase integration" trick to represent 1204(1 as an integral of some determinant, but there will be two 
phase variables, since there are 3 types of inequivalent brackets: 

S b (Ki,...,KM;d,r) = (r-l)!(M-r)!(-l) d - 1 x 

5 b (V',0, Kl,«l) g b (lp,(l),K 1 ,K2) ... g b (ip,(f>,K 1 ,K M ) 
g b (ip,<p,K 2 ,K 1 ) g b (ip,<f>, K 2 ,K 2 ) ... g b (lp,(j>,K2,KM) 



(r-l>0# f n d^ p -^ 



2tt 



2vr 



det 



g b (lp,(f),KM,Kl) g b {^, </>, K M ,K 2 ) ... g b (i),(f),KM,K M ) 



where 



5 b (</>, 0, j, = e^(d, j, I) + e^gl(d,j, I) + g 2 b (d,j, I). 



(206) 



(207) 



After integration over (j), determinant in 1(206(1 has terms up to e % ^ M therefore integration over ip and summation 
according to l(201|) are equivalent to substitution ip$ = : 



S b ( Ku ...,K M ;d) = (N - M)\(M - l)l(-l) 



d-l 



2tt 



2tt 



-e x 



det 



ff 6 (0, <j>, Kl,Kl) g b (0,(f>,Kx,K2) ... g b (0,<f>,Ki,K M ) 
3 b (0, 0, «2,«i) # b (0, 0, k 2 ,k 2 ) ■■■ 5 b (0, 0, k 2 ,k m ) 

3 b (0, 0, Kmj«i) 3 b (0, 0, k m ,k 2 ) ... g b (0, 0, km, k m ) 



(208) 



Integral over can be simplified further, since the determinant in ((208|l has a form Ao+vlie 1 ^. The form above follows 
from the fact that a part of the matrix which depends on e*^ has a rank 1 and the formula for the determinant of 
the sum of the matrices(see page 221 o&>). Let's for a moment introduce a notation z = e 1 ^ . Integration over with 
a weigh e~ lip extracts the term A%, which can be alternatively written as a difference between two determinants, one 
when z = 1 and the other when z — (gf(ip, j, is given by 1(192(1 ) 



S"(/si, k m ; d) = (N - M)!(M - l^-lf" 1 det 



(JV-M)!(M- l)!(-l) d_1 det 



g b (0,0,Ki,Ki) gr b (0,0, ki,k 2 ) ... g b (0, 0,ki,k m ) 

5 b (0,0,K 2 ,K!) ff b (0,0,K 2 ,K 2 ) ... 5 b (0,0,K 2 ,K M ) 

3 b (0, 0, k m , ki) g b (0, 0, k m , k 2 ) ... g b (0, 0, k m , k m ) _ 

5^(0, ki,ki) g / (0, ki,k 2 ) ... 3 / (0,ki,k m ) 

g / (0, k 2 ,ki) g / (0, k 2 ,k 2 ) ... g f (0,K 2 ,K M ) 



g f (0,K M ,Ki) g / (0, k m ,k 2 ) ... 5^(0, k m , k m ) 



.(209) 



We note, that a similar trick is explained on the page 609 ofS. After that, summation over different Ki can be 
performed similarly to the case of fermions: 



S b (d;p;T) = (N - M)\{M - 1)!(-1) 



e i(N-M)6 x 
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det 



e" + /(l) 5 6 (0, 0,1,1) /(iy(0,0,l,2) 
/(2). 9 h (0, 0,2,1) e' 9 + /(2) 5 b (0,0,2,2) 



f(N)g b (0,0,N,l) f(N)g b (0,0,N,2) 
where S f (d;p;T) is defined in (|196jl . 



/(l) 5 b (0,0,l,7V) 
/(2) 3 b (0,0,2,7V) 

e' 9 + /(iV)5 b (0,0,iV, AT) 



N -M 
M 



Sf(d;p;T), (210) 



VII. EXPERIMENTAL CONSIDERATIONS AND CONCLUSIONS 



In this section we will consider in detail possible ways to realize the system under investigation in experiments with 
cold atoms. 

An array of one dimensional tubes of cold atoms along x direction has been realized experimentally using strong 
optical lattices in two dimensionsi&Si^S^ y and z. The large number of tubes provides a good imaging quality, 
but the number of atoms and the ratio between bose and fermi particle numbers varies from tube to tube, and may 
complicate the interpretation of the experiments (one of the ways to fix the ratio between bose and fermi numbers 
for all tubes will be discussed later). In addition, due to harmonic confinement along the axis of the tube, bose 
and fermi densities vary within each tube, which causes non-homogeneous broadening of the momentum distribution. 
Alternatively, single copies of one dimensional mixtures with constant densities along the axis can be realized in micro 
traps on a chip£», or using cold atoms in a Id box potential^. Here we will mostly concentrate on a realization of Id 
system using strong 2D optical lattice in y and z directions. 

First of the conditions mj = to/, is approximately satisfied for isotopes of the atoms, and one can expect our 
theory to be valid with high accuracy for them. Some of the promising candidates are 39 ( 41 ^K— 40 K—, 171 Yb+ 172 Yf£L, 
and 86 ( 84 )i?6 — 87 ( 85 ) Rb&L. Different isotopes of potassium have already been cooled to quantum degeneracyiSi by 
sympathetic cooling with Rb. There is another way to satisfy the first condition of J5J) using already available degenerate 
mixtures^. If one uses an additional optical lattice along the x direction with filling factors much smaller than one, then 
is an effective Hamiltonian describing this system with the effective masses determined by the tunneling, similarly 
to a recent realization of Tonks- Girardeau gas for bosons 8 . Finally, we note that one can realize experimentally the 
model, which has the same energy eigenvalues as using a mixture of two bosonic atoms (see next paragraph). If 
one chooses two magnetic sub-levels of the same atom, equality of masses will be satisfied automatically. 

Second of the conditions J2J, gbb = 9bf > 0, can also be satisfied in current experiments, using a combination of sev- 
eral approaches. First, one can use Feshbach resonances to control the interactions: this is particularly straightforward 
for Li — Na of K — Rb mixtures, where resonances have already been observed experimentally^*. Second, we point out 
that it is sufficient to have equal (positive) signs for the two scattering lengths, but not necessarily their magnitudes. 
Well away from confinement induced resonances^, ID interactions are given by gbb = 2htOb±o-bb, 9bf = ^^^b±^f±dbf, 
where ujb±, Wf± are radial confinement frequencies, and abb, a-bf are 3D scattering lengths. For a fixed value of dbb/a>bfi 
one can always choose the detuning of the optical lattice laser frequencies in such a way that gbb = 9b f- After that, one 
can vary the intensity of the y, z optical lattice beams and change g, while always being on the integrable line of the 
phase diagram. Combination of these two approaches to control ID interactions gives a lot of freedom for experimental 
realization of equal one dimensional interactions. Finally, lets describe how to realize the bosonic model, which has 
the same eigenvalues as the model l|T)l. Bosonic system is characterized by 3 interaction parameters, gn, 322, 312- If one 
tunes gn to +00, then bosons of type 1 get " fermionized" within the same type, and the model will be equivalent in 
terms of energy spectrum, density profiles and collective modes to Note, however, that single-particle correlation 
functions will be different, and the results of sections IVl and fvTI (except for IVI B|) are not applicable. This general 
equivalence between bose-bose and bose- fermi models is valid for any ratio between (722 and 512. One can push this 
result even further, by tuning 322 to +00. In this case eigenstates of Q are equivalent to spin— 1/2 fermi systemSSiSiiSI, 
and some predictions for those systems can be applied for bosons. 

Detection of the properties of the system may be hindered by the fact, that both number of atoms and relative 
fraction of bosons a vary from tube to tube. However, one can use Feshbach resonances to fix the boson fraction to be 
a = 1/2 in each tubaSi. To do this, one can use Feshbach resonance for bose- fermi scattering to adiabatically create 
molecules before loading the mixture in strong y, z optical lattice. If one gets rid of unpaired atoms at this stage, 
switches on y, z optical lattice, and adiabatically dissociates the molecules, boson fraction will be fixed in each tube 
to be a = 1/2. Most of our figures have been calculated for this particular boson fraction. Our results in harmonic 
traps are presented as functions of 70 = ragj (h 2 no), where no is a total density in the center of a one dimensional 
trap, and 70 is the Lieb-Liniger— parameter in the center of the trap. 70 3> 1 corresponds to a strongly interacting 
regime, rig varies from tube to tube, and to be able to compare theoretical predictions precisely with experiments, 
one should be able to have an optical access to regions where variation of no is small. 
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Most of our experimental predictions, except for those in section IVT1 deal with zero temperature case. Experimen- 
tally, one needs to verify the quantum degeneracy of the gases in ID regime. A possible way to identify the onset of 
quantum degeneracy is based on density profiles^. In Figs. |2and|21we show the density profiles at zero temperature 
for weak and strong interactions, when the harmonic confinement frequency loq is the same for bosons and fermions. 
In both cases, only central part is occupied by bosons, and outer shells consist of fermions only. In addition, for 
the strong interactions fermi density develops a strong peak at the edge of bosonic cloud. When the interactions 
are not strong (70 < 1), one can estimate the temperature at which quantum effects become important for ground 
state density profile to be of the order of Nfkoo, where N is the total number of atoms in a tube. In the strongly 
interacting regime (70 3> 1), however, situation is very different. There are two temperature scales in the problem: 
Ef = (nhn ) 2 1 (2m) , and -Ej/70 *C E®. As the temperature goes up from to ~ Ej/^ , density profile changes as 
shown in figure 1141 and the peak in the fermion density disappears. However, total density profile doesn't change 
much as long as T <C E®. This effect can be qualitatively understood as the demonstration of the " fermionization" of 
the bose-fermi cloud, as will be explained in the next paragraph. 

First, lets consider the case without a harmonic potential. When interactions are strong, bosons tend to avoid 
fermions and other bosons. Whenever coordinates of any two particles coincide, wavefunction is close to 0. Effectively, 
the gas is mutually " fermionized" , and the ground state energy of the system is close to the ground state energy of 
the pure noninteracting fermi gas with a density equal to the total density of bosons and fermions. Dependence of the 
energy on the relative density (or boson fraction a) appears only in the next order in 1/7 expansion, and two first 
terms in this expansion are given by l|54|l . Since dependence of the energy on boson fraction a is 7 3> 1 times smaller 
than dependence on total density, the "quantum degeneracy" temperature for relative density excitations is also 7 
times smaller than quantum degeneracy temperature for fermions with density n, hence it is ~ Ef/'y. When harmonic 
trap is present at T = 0, relative density distributes itself to minimize the total energy. As temperature becomes of 
the order of several E^/jq, almost all relative density modes get excited, and boson fraction becomes uniform along 
the trap. Total density modes are still not excited, since their quantum degeneracy temperature is E®, and therefore 
the total density profile doesn't change much. Temperature Ef/jo, is important not only for density distribution, but 
also for correlation functions, as will be discussed later. 

Knowledge of the exact dependence of the energy as the function of densities and interactions allows to investigate 
not only the static properties, but also dynamic behavior. In section ItVl we developed a two-fluid hydrodynamic ap- 
proach to calculate the frequencies of collective oscillations. In the strongly interacting limit we predict the appearance 
of low-lying modes, with a frequency scaling as ~ ^o/a/To- These modes correspond to "out of phase" oscillations 
of bose and fermi clouds that keep the total density approximately constant. These modes can be understood as 
follows: due to fermionization effects discussed in previous paragraph, for 70 3> 1 the energetic penalty for changing 
the relative density of bosons and fermions is small, and hence it doesn't cost too much energy to create "out of 
phase" oscillations that don't change the total density. Dependence of the frequencies of low-lying oscillations with 
small quantum numbers on overall boson fraction in a tube is shown in figure |H1 In addition to low lying "out of 
phase" oscillations, the cloud has "in phase" oscillations, with the frequencies w„ = hloq, similarly to Tonks-Girardeau 
gas of bosons 3 ?. These modes have frequencies considerably higher than "out of phase" modes, and are not shown 
in figure |H1 One can excite any of these excitations by adding a perturbation of the matching frequency, similarly to 
what has been done to bosons in2£. A different manifestation of the slow "out of phase" dynamics can be observed 
looking at the evolution of density perturbations: initial perturbation will split into fast "in phase" part, moving 
at fermi velocity, and slow "out of phase" part. This is similar to "spin-charge separation", proposed for fermiSi or 
bose^S spin 1/2 mixtures. When interactions are not strong (70 < 1), one can obtain frequencies of all modes using 
mean- field energy. Figure^ shows the dependence of frequencies for equal number of bosons and fermions (a = 1/2) 
on 7o. Even in mean field regime frequency of "out of phase" oscillations gets smaller as interactions get stronger. 
Already for 70 « 1 results for 70 3> 1 extrapolate mean-field results very well. 

Finally, lets discuss theoretically the most interesting and sensitive measure of the correlations, single particle 
correlation function, considered in sections and IVII Fourier transform of the single particle correlation func- 
tion is an occupation number, and it can be measured experimentally using Bragg spectroscopy^ or time of flight 
measurements^. We can calculate these correlation functions in strongly interacting regime under periodic boundary 
conditions for any temperatures. At zero temperature bose momentum distribution has a sing ularity ltToT|) at k = 
reminiscent of BEC in higher dimensions, and its strength is controlled by Luttinger liquid parameter Kb, which 
depends only on boson fraction for strong interactions. For fermions, momentum distribution has a lot of interesting 
features. At zero temperature, several momentum distributions are presented in figs, ^fl El El O ne sees, that due 
to strong interactions, fermi step at kf gets smeared out even at T — 0, and n'[k) is considerably different from 
at wave vectors far away from kf. However, total change of n/(fc) as one crosses kf is quite large. In addition, n/(fc) 
develops an extra singularity^ at kf + 2kb, and the strength of this singularity is higher for small boson fractions. As 
the temperature rises, momentum distribution changes considerably in the region of low temperatures of the order 
of Ef/j, and its evolution as a function of temperature is shown in figure IT51 For Ef <C T <C Ef, one enters so 
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r 2tc^ 

FIG. 16: Momentum distribution for fermions after averaging by inhomogeneous density profile for harmonic confinement. 
Results are shown for T = 0, (thick line), E^/jo <C T <C E°j (normal line) and for the same number of noninteracting fermions 
(dashed). Overall number of bosons in a trap equals the total number of fermions, no is the total density in the center of the 
trap, E® = (Trhno) 2 / (2m) , 70 2> 1 is the Lieb-Liniger parameter in the center of the trap, and 2xf is the total size of the cloud. 
In the range of the temperatures ~ E® /70 fermi correlation function changes considerably due to transition from true ground 
state to "spin disordered" regime. 

called "spin disordered" regim e 46 ' 4 ? , where singularity at kf gets completely washed out, and for equal densities of 
bosons and fermions momentum distribution gets almost twice as wide compared to T <C E*j /j. A strong change of 
the momentum distribution in a small range of temperatures can be used to perform a thermometry at very small 
temperatures. To verify experimentally exact numerical correlation functions one needs to work with systems at con- 
stant densities along x direction. Such constant density can be achieved in experiments with micro traps^, or in 2D 
arrays of tubes, if one makes a very shallow harmonic confinement, and creates strong box-like impenetrable potential 
at the sides of the tubes with the help of additional lasers. If the system is in harmonic trap, lots of the features of 
correlations themselves (i.e. singularity at kf + 2kf>) get washed out due to averaging over inhomogeneous density 
profile 63 . However, the averaged correlation function still shows significant change in the region of temperatures of the 
order of and the results for T <C £7/70 and Ef /70 <C T <C are shown in fig. ^] The point where (k) 

has a discontinuous derivative for T — corresponds to the fermi wavevector for the maximal density of fermions (at 
the edge of the bosonic cloud). For comparison, we also show N*(k) for the same number of fermions in the same 
trap for noninteracting case. 

In conclusion, we presented a model for interacting bose-fermi mixture in ID, which is exactly solvable by Bethe 
ansatz technique. We obtained the energy numerically in the thermodynamic limit, and used it to prove the absence 
of the demixing under conditions J2J, contrary to prediction of a mean field approximation. Combining exact solution 
with local density approximation (LDA) in a harmonic trap, we calculated the density profiles and frequencies of 
collective modes in various limits. In the strongly interacting regime, we predicted the appearance of low-lying 
collective oscillations which correspond to the counterflow of the two species. In the strongly interacting regime we used 
exact wavefunction to calculate the single particle correlation functions for bosons and fermions at zero temperature 
under periodic boundary conditions. We derived an analytical formula, which allows to calculate correlation functions 
at all distances numerically for a polynomial time in system size. We investigated numerically two strong singularities 
of the momentum distribution for fermions at kf and kf + 2kb- We extended the results for correlation functions 
for low temperatures, and calculated correlation functions in the crossover regime from T = to "spin disordered" 
regime. We also calculated the evolution of the density profile in a harmonic trap at small nonzero temperatures. We 
showed, that in strongly interacting regime correlation functions change dramatically as temperature changes from 
to a small temperature ~ £7/7 <§C Ef, where Ef = (irfin) 2 / (2m) , n is the total density and 7 is the Lieb-Liniger 
parameter. Finally, we analyzed the experimental situation, proposed several ways to implement the exactly solvable 
hamiltonian and combined the results for correlation functions with LDA. 

We thank M. Lukin, L. Mathey, G. Shlyapnikov, D.Petrov, P.Wiegmann, C. Menotti and DW. Wang for useful 
discussions. This work was partially supported by the NSF grant DMR-0132874. 




35 



APPENDIX A 

In this appendix we will prove that all solutions of equations (|28[1 - <|29[) 



N 



n kj~A a + ic/2 
ib i -A Q -ic/2 =1 ' a = {1 '-' M} ' 

^'- A ^ ic/2 ,i = {l,..^} 



M 

n 



_ kj — Ap — ic/2 



are always real. This is a major simplification for the analysis of the excited states compared to spin- 
systems, where one has to consider complex solutions^. 

Suppose that solutions of (|A1|) - (|A2|) are complex numbers, such that 

inf Imfcj = k~ < suplmfcj = k + , 
inf Im A Q — A~ < sup Im A a = A + . 



We need to prove that k 
First, lets prove that 



k+ = A- = A+ = 0. 



Suppose that (|A5|1 is not valid, i. e. 



Then 



k~ < A", 
A+ < k+. 



3 a : Im kj — Im A Q > V j. 

kj — A a + ic/2 



kj — A a + ic/2 



>lVj, 



and absolute value of the lhs of eq. IjAlfl is bigger than 1, which contradicts the equation. Equation 
proven similarly. 

Now, lets prove that 

k+ < 0, 
k" > 0. 

These equations together with ljA5(l - l|A6(l would imply k~ = k + = A - = A + = 0. 
Suppose that ljA9() is not valid, i. e. 3 j : Im kj = k + > 0. From (|A6|I it follows that 



therefore 



Imkj - ImAp > V (3, 
kj - Ap + ic/2 



kj - Ap + ic/2 



> 1 Vj9, 



and absolute value of the rhs of equation I)A2J| is not smaller than 1. On the other hand, by assumption 
equation is smaller than 1 : 



|e*^l = e~ k+L < 1. 



(Al) 

(A2) 
fcrmion 



(A3) 
(A4) 



(A5) 
(A6) 



(A7) 

(A8) 
6l can be 



(A9) 
(A10) 



(All) 

(A12) 
lhs of this 

(A13) 



Contradiction proves the validity of (|A9(I . and l|A10(l can be proven similarly. 
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